load packages:

library(googlesheets4) #scrapes google sheets
library(dplyr) # dataframe manipulation
library(ggplot2) # plotting package
library(lubridate)
library(tidyr)
library(readr)
library(scales)
library(growthrates) 
library(patchwork)
library(grid)
library(qiime2R)
library(phyloseq)
library(viridis)
library(magrittr)
library(vegan)
library(stringr)
library(ape)
library(RColorBrewer)

Create a color pallete

pp_palette <- c(
  "#7D8F69", "cornflowerblue", "#F5C45E",  "plum",
  "hotpink4", "lightsalmon", "steelblue", "indianred", "springgreen4", 
   "#102E50", "maroon", "lightgreen", "lightblue", "gold",
   "#FF0099", "#A26769", "blue", "red", "lightpink", 
  "#71816D", "#FFE6E6", "#FF66CC", "#A3B18A", "salmon", "magenta","orchid","deepskyblue", "tomato",
    "cyan", "hotpink", "yellow", "orange", "navy", "darkred",
   "#C0A98E", "#7A9E9F", "#A1869E", "#BFD8B8", "#F7DD72","green", "darkgreen","purple", "coral", "#FFA500", 
  "#FFE5AD",
  "#CDB4DB", "pink","dodgerblue1","green2", "brown4", "darkviolet", "azure2", "aquamarine"
)
scales::show_col(pp_palette)

1 Characterization of bacterial communities in P. bahamense cultures

1.1 Import sequence data from QIIME 2:

Amplicons were processed on a QIIME 2 pipeline, which includes DADA denoising. Seqeunces were classified with Native Bayers Classifier

ps<-qza_to_phyloseq(features="table.qza")

1.2 Import and Format Sample Metadata:

meta <- read.csv("SampleMetaData.csv")
meta$SampleID <- paste0(meta$SampleID, "_6283")
row.names(meta)<- meta$SampleID
meta$TreatRep <- paste(meta$PCRDate_Treat, meta$SampleInfo_Rep,  sep = "_")
#keep only the P. bahamense project
pyro_only <- meta[grepl("^Pyro", meta$Project, ignore.case = TRUE), ]
META <- sample_data(pyro_only)

1.3 Import and Format Taxonomy Data:

taxonomy <- read.csv("taxonomy.csv", stringsAsFactors = FALSE)

names(taxonomy) <- c("row", "tax", "Confidence") 
row.names(taxonomy) <-taxonomy[[1]] #move the feature ID column to become the row names
taxonomy <- taxonomy[,(-1)] #delete the feature ID  column 
taxonomy <-  separate(taxonomy, tax, c("Domain","Phylum", "Class", "Order", "Family", "Genus", "Species", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14"), sep = ";", fill = "right")
taxonomy <- taxonomy[,c(1:7)] #keep first seven levels

taxonomy$D0 <- with(taxonomy, ifelse(Order == " o__Chloroplast", "Chloroplast", "Bacteria")) #create column for ordering

# Reorder columns and turn into phyloseq tax_table
col_order<- c("D0", "Domain","Phylum", "Class", "Order", "Family", "Genus", "Species")
taxonomy<- taxonomy[, col_order]
taxmat <- as.matrix(taxonomy)
TAX = tax_table(taxmat)

1.4 Build a phyloseq object:

#Keep only Pyro samples
ps <- prune_samples(rownames(pyro_only), ps)

#Remove OTUs that have zero counts in all Pyro samples
ps <- prune_taxa(taxa_sums(ps) > 0, ps)
ps <- merge_phyloseq(ps, TAX, META)

1.5 Domain-Level Relative Abundance Bar Plot:

#Keep only Bacteria + Chloroplast
ps<- subset_taxa(ps, D0 %in% c("Bacteria", "Chloroplast") )

#Convert counts to relative abundance (%)
psra<- transform_sample_counts(ps, function(OTU) 100* OTU/sum(OTU)) 

#Group OTU by DO level
glomD1<- tax_glom(psra, 'D0')

#plot RA
taxabarplot <- plot_bar(glomD1, x = "SampleID", fill = "D0") +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    text = element_text(size = 14),
    axis.text.x = element_text(angle = 90)
  ) +
  ylab("Relative Abundance (%)") +
  xlab("") +
  scale_fill_manual(values = c("lightgrey", "#7BB03B")) +
  facet_grid(~Project, scales = "free_x")


taxabarplot

1.6 Remove chloroplast sequences from data:

ps_nochloro <- subset_taxa(ps, Order != " o__Chloroplast" & Domain != "Unassigned")

1.7 Calculate observed richness:

plugin <- ps_nochloro %>%
            estimate_richness(measures = "Observed") %$% Observed

#Pull project metadata field from phyloseq object
Project <- ps_nochloro %>% sample_data %$% Project

#Combine richness + project into one dataframe
richness<- data.frame(plugin, Project)
names(richness) <- c("richness", "Project")


richness %>%group_by(Project) %>% summarize(mean = mean(richness), min = min(richness), max = max(richness))

1.8 Relative abundance and sums the counts of all OTUs that belong to the same Order:

ps_nochloro_RA<- transform_sample_counts(ps_nochloro, function(OTU) 100* OTU/sum(OTU))

ps_nochloro_RA_glomO<- tax_glom(ps_nochloro_RA, 'Order')

1.9 Bar plot at the order level:

# subset by pyro samples
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")

# Remove low abundance orders
pyro = filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
                       
taxabarplot <- plot_bar(pyro, x= "SampleID", fill = "Order") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Order), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

1.10 Bar plot at family level:

ps_nochloro_RA<- transform_sample_counts(ps_nochloro, function(OTU) 100* OTU/sum(OTU))

#Group OTU by family and sum abundances
ps_nochloro_RA_glomF<- tax_glom(ps_nochloro_RA, 'Family')

#Keep only the Pyro samples
pyro <- subset_samples(ps_nochloro_RA_glomF, Project == "Pyro")

#Remove low abundance orders
pyro = filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
                       
taxabarplot<-plot_bar(pyro, x= "SampleID", fill = "Family") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Family), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

1.11 Filter Samples and Generate Family-Level Bar Plot:

#Filter SampleID
pyroE <- subset_samples(pyro, !(SampleID %in% c("L14_6283", "L15_6283", "L17_6283", "L18_6283")))

#Remove low abundance orders
pyroE = filter_taxa(pyroE, function(x) sum(x) > 1, TRUE)
                       
taxabarplot<-plot_bar(pyroE, x= "SampleID", fill = "Family") +  scale_y_continuous(expand = c(0, 0)) + ggtitle("")  + theme(legend.title=element_blank()) + geom_bar(aes(fill=Family), stat="identity", position="stack", width =0.9) +theme_classic() + theme(text = element_text(size=14))+theme(axis.text.x = element_text(angle = 90)) + xlab("Sample") + ylab("Relative Abundance (%)") + theme(text = element_text(size=14)) + scale_fill_manual(values= pp_palette)+ xlab("")

taxabarplot

1.12 Determine Unique OTU’s and abundances in P. bahamense data:

pE_m <- psmelt(pyroE)

unique(pE_m$OTU)
## [1] "12a986d7ed4f3be699daa0feaba41d8a" "07ae9a3ebf44d4a821d001debe7d9b0c"
## [3] "e9034a72595933ebd5c5df08eb4a08fc"
  1. 12a986d7ed4f3be699daa0feaba41d8a

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

12a986d7ed4f3be699daa0feaba41d8a GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

uncultured sanger sequence from Mangrove Sediment in Hong Kong - GenBank: MH091208.1, MH091199.1 (2 best hits, 100% ) Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs

Diversity of salt marsh prokaryotes, M.A. Moran, uncultured Hyphomicrobiaceae bacterium - lon=81.2797W, lat=31.3884N; sediment 14-16cm collected on Feb 01, 2002, Sapelo Island Microbial Observatory Dean Creek Marsh sampling site” (3rd best hit 99%)

pE_m %>%  filter(OTU == "12a986d7ed4f3be699daa0feaba41d8a") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
  1. e0169d54945a8a584037e96365ec7bb3

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

e0169d54945a8a584037e96365ec7bb3 GCTGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria - Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - isolation_source=“mangrove sediments” in Hong Kong, 99%

pE_m %>%  filter(OTU == "e0169d54945a8a584037e96365ec7bb3") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
  1. aa0c627da86b1ddbba980c9dfd6ac380

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

aa0c627da86b1ddbba980c9dfd6ac380 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGACTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria - Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - Hong Kong isolation_source=“mangrove sediments” , 99%

pE_m %>%  filter(OTU == "aa0c627da86b1ddbba980c9dfd6ac380") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
  1. f6a0f1f2d2a083164d1fb145b000b96b

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

f6a0f1f2d2a083164d1fb145b000b96b GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGTGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGAGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 Uncultured bacteria Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs - isolation_source=“mangrove sediments” - 99%

pE_m %>%  filter(OTU == "f6a0f1f2d2a083164d1fb145b000b96b") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
  1. ff08ec2d8a67e142bc459567275bb888

o__Rhizobiales f__Hyphomicrobiaceae g__Filomicrobium s__uncultured_Alphaproteobacteria

ff08ec2d8a67e142bc459567275bb888 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAACTCTTTTGGCGGGGACGATAATGACGGTACCCGCAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGACTGGTCAGTTGGGGGTGAAATCCCAGGGCTCAACCCTGGAACTGCCTCCAATACTGCCAGTCTTGAGTCCGAGAGAGGTGACTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

MH091208.1 uncultured bacteria Diversity and dynamics of microbial community structure in different mangrove, marine and freshwater sediments during anaerobic debromination of PBDEs Hong Kong isolation_source=“mangrove sediments” 99%

EU488009.1 uncultured bacteria Characterization of the lucinid bivalve-bacteria symbiotic system: the significance of the geochemical habitat on bacterial symbiont diversity and phylogeny isolation_source=“siliciclastic sedment from Thalassia sea grass bed” 99%

other marine sediments, estuary sediments

pE_m %>%  filter(OTU == "ff08ec2d8a67e142bc459567275bb888") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))

  1. 07ae9a3ebf44d4a821d001debe7d9b0c

c__Acidimicrobiia o__Microtrichales f__uncultured g__uncultured NA

07ae9a3ebf44d4a821d001debe7d9b0c GCAGCAGTGGGGAATCTTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGCGGGAAGACGGCCTTCGGGTTGTAAACCGCTTTCAGCAGGGACGAAATTGACGGTACCTGCAGAAGAAGCTCCGGCCAACTACGTGCCAGCAGCCGCGGTAAGACGTAGGGGGCGAGCGTTGTCCGGAATCATTGGGCGTAAAGGGCTCCTAGGTGGTTCAGTAAGTCGACTGTGAAAATCCAAGGCTCAACCTTGGGACGCCAGTCGATACTGCTGTGACTCGAGTTCGGTAGAGGAGTGTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCAACGGCGAAGGCAGCACTCTGGGCCGATACTGACACTGAAGAGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

KM840989.1 uncultured bacteria - Microbial diversity of indigenous bacteria in a 129I contaminated groundwater plume at the Hanford Site, Washington - Sanger, isolation_source=“iodine (I-129) contaminated groundwater” 100%

HF558551.1 uncultured bacteria - Iron- and Sulphur- cycling bacteria mobilize copper in a multiple extreme mine tailings in the Atacama Desert, Chile - isolation_source=“tailing material” - 100%

JQ427269.1 uncultured bacteria - Bacterial diversity in an alkaline saline soil spiked with anthracene, samger sequenced, isolation_source=“soil”, 100%

JQ665390.1 uncultured actinobacterium, Diversity of unculturable Actinomycetes in coastal wetlands of the Yellow River estuary, isolation_source=“soil sample from coastal wetlands”, 99%

pE_m %>%  filter(OTU == "07ae9a3ebf44d4a821d001debe7d9b0c") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))
  1. 9b8accad19c30d7a69a6290553111166

c__Acidimicrobiia o__Microtrichales f__uncultured g__uncultured NA

9b8accad19c30d7a69a6290553111166 GCTGCAGTGGGGAATCTTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGCGGGAAGACGGCCTTCGGGTTGTAAACCGCTTTCAGCAGGGACGAAATTGACGGTACCTGCAGAAGAAGCTCCGGCCAACTACGTGCCAGCAGCCGCGGTAAGACGTAGGGGGCGAGCGTTGTCCGGAATCATTGGGCGTAAAGGGCTCCTAGGTGGTTCAGTAAGTCGACTGTGAAAATCCAAGGCTCAACCTTGGGACGCCAGTCGATACTGCTGTGACTCGAGTTCGGTAGAGGAGTGTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCAACGGCGAAGGCAGCACTCTGGGCCGATACTGACACTGAAGAGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

KM840989.1 uncultured bacteria - Microbial diversity of indigenous bacteria in a 129I contaminated groundwater plume at the Hanford Site, Washington - isolation_source=“iodine (I-129) contaminated groundwater” - 99%

pE_m %>%  filter(OTU == "9b8accad19c30d7a69a6290553111166") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))

  1. e9034a72595933ebd5c5df08eb4a08fc

o__Oceanospirillales f__Alcanivoracaceae1 g__Alcanivorax s__Alcanivorax_venustensis

e9034a72595933ebd5c5df08eb4a08fc GCAGCAGTGGGGAATCTTGGACAATGGGGGCAACCCTGATCCAGCCATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTCAGCAGGGAGGAAGGCTTACCCCTAATACGGGTGAGTACTTGACGTTACCTGCAGAAGAAGCACCGGCTAATTTCGTGCCAGCAGCCGCGGTAATACGAAAGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGCGGTGTGTTAAGTCGGATGTGAAAGCCCAGGGCTCAACCTTGGAATTGCATCCGATACTGGCACGCTAGAGTGCAGTAGAGGGAGGTGGAATTTCCGGTGTAGCGGTGAAATGCGTAGAGATCGGAAGGAACACCAGTGGCGAAGGCGGCCTCCTGGACTGACACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

PP516259.1 organism=“Alloalcanivorax venustensis” Exploration and conservation of bacterial community from the Arabian Sea seamount isolation_source=“Arabian Sea water”, 100%

pE_m %>%  filter(OTU == "e9034a72595933ebd5c5df08eb4a08fc") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))

  1. c02844d75d512c1510de7334b6687731

o__Rhizobiales f__Hyphomicrobiaceae g__uncultured s__uncultured_bacterium

c02844d75d512c1510de7334b6687731 GCAGCAGTGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCCATGCCGCGTGAGTGACGAAGGCCTTAGGGTTGTAAAGCTCTTTTGGCGGGGAAGATAATGACGGTACCCGCAGAATAAGCTCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGAGCTAGCGTTGTTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGATTTGTTAGTCAGGGGTGAAATCCCGGGGCTCAACCCCGGAACTGCCTTTGATACTGCAAATCTCGAGTCCGAGAGAGGTGGGTGGAATTCCTAGTGTAGAGGTGAAATTCGTAGATATTAGGAAGAACACCGGTGGCGAAGGCGGCCCACTGGCTCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATAC

OK235762.1 Uncultured Pedomicrobium sp. isolation_source=“date palm rhizosphere soil” Saudi Arabia - 99%

HQ697761.1 Uncultured bacteria isolation_source=“hydrocarbon contaminated saline-alkali soil” China - 99%

99% to several other contaminated soil sites (real and experimental)

KP098952.1 Uncultured bacteria - The microbiome of methanol-utilizing denitrification systems contains new bacterial groups - isolation_source=“denitrification bioreactor” - 99%

pE_m %>%  filter(OTU == "c02844d75d512c1510de7334b6687731") %>%  summarize(mean = mean(Abundance), min = min(Abundance), max = max(Abundance))

1.13 Plot shannon diversity and relative abundance at order level:

# -------------------------------
# Subset Pyro samples and filter
# -------------------------------
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")
pyro <- subset_samples(pyro, Treatment != "AB")
pyro <- filter_taxa(pyro, function(x) sum(x) > 1, TRUE)
pyro <- subset_samples(pyro, 
                       Treatment %in% c("1003", "Replete") & 
                       TreatRep %in% c("1003_1", "1_2"))

# Set factor levels for Notes_Filter
pyro@sam_data$Notes_Filter <- factor(as.character(pyro@sam_data$Notes_Filter),
                                     levels = c("10", "2"))

# Clean Order names in tax_table
tax_table(pyro)[, "Order"] <- gsub(" o__", "", tax_table(pyro)[, "Order"])

# -------------------------------
# Plot relative abundance barplot
# -------------------------------
barplot <- plot_bar(pyro, x = "TreatRep", fill = "Order") +
  geom_bar(aes(fill = Order), stat = "identity", position = "stack", width = 1) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_manual(values = pp_palette) +
  facet_grid(Notes_Filter ~ Treatment, scales = "free") +
  scale_x_discrete(
    expand = c(0.5, 0.05),
    labels = c(
      "1003_1" = "IRL",
      "1_2"    = "OTB"
    )
  ) +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    legend.position = "right",
    legend.justification = "center",
    legend.box = "horizontal",
    axis.text.x = element_text(angle = 0, vjust = 0.5, size = 16, color = "black"),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.title.x = element_text(size = 20, color = "black"),
    axis.title.y = element_text(size = 20, color = "black"),
    axis.ticks.x = element_line(color = "black", linewidth = 0.9),
    axis.ticks.y = element_line(color = "black", linewidth = 0.9),
    strip.text.x = element_blank(),  # remove top facet labels
    strip.text.y = element_text(size = 20, color = "black"), # row facets
    text = element_text(size = 20, color = "black")
  ) +
  ylab(expression("Relative Abundance (%)")) +
  xlab("")

barplot

# -------------------------------
# Calculate Shannon index
# -------------------------------
pyro <- subset_samples(ps_nochloro_RA_glomO, Project == "Pyro")
pyro <- subset_samples(pyro, TreatRep %in% c("1003_1", "1_2", "1_1", "1_3"))
pyro <- subset_samples(pyro, Treatment != "AB")

shannon_df <- estimate_richness(pyro, measures = "Shannon")

# Extract metadata and combine with diversity
sample_meta <- as(sample_data(pyro), "data.frame")
shannon_df <- cbind(shannon_df, sample_meta)

# -------------------------------
# Plot Shannon diversity
# -------------------------------
shannon_plot <- ggplot(shannon_df, aes(x = Treatment, y = Shannon, color = Notes_Filter)) +
  
  # Color mapping
  scale_color_manual(values = c("10" = "navy", "2" = "maroon")) +
  
  # Points
  geom_point(size = 5) +
  
  # Custom x-axis labels
  scale_x_discrete(labels = c(
    "1003" = "IRL",
    "Replete" = "OTB"
  )) +
  
  # Classic theme
  theme_classic() +
  
  # Axis labels
  ylab("Shannon Diversity Index") +
  xlab("") +
  labs(color = "Filter Size") +  

  # Theme adjustments
  theme(
    legend.title = element_text(size = 20),
    axis.text.x = element_text(angle = 0, vjust = 0.5, size = 16, color = "black"),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.title.x = element_text(size = 20, color = "black"),
    axis.title.y = element_text(size = 20, color = "black"),
    axis.ticks.x = element_line(color = "black", linewidth = 0.9),
    axis.ticks.y = element_line(color = "black", linewidth = 0.9),
    strip.text.x = element_blank(),  # remove top facet labels
    strip.text.y = element_text(size = 20, color = "black"), # row facets
    text = element_text(size = 20, color = "black"))

# Print Shannon plot
shannon_plot

1.14 Figure 2

bo <- shannon_plot + barplot
bo

Fig. 2 Alpha diversity and composition of bacterial communities in P. bahamense cultures. A Shannon diversity index for particle-associated (>10.0-µm size fraction; blue) and free-living (0.2–10.0-µm size fraction; red) bacterial communities in cultures of FWC1003 and FWC1006 P. bahamense isolated from Indian River Lagoon and Old Tampa Bay, respectively. B Relative abundance of bacterial orders in size-fractioned communities in P. bahamense cultures.

2 Evaluating the B-vitamin requirements of P. bahamense

Chlorophyll-a fluorescence (RFU) data from P. bahamense cultures were used to visualize growth patterns, quantify maximum biomass, and estimate specific growth rates under media conditions deficient in media lacking one targeted B vitamin. Specific growth rates were computed using the “alleasy_linear” function of the “growthrates” package in R Studio (Hall et al., 2014; RStudio Team, 2020).

2.1 Import the Dataset:

fluor<-as.data.frame(read_csv("fluor.csv"))
fluor$Time <- as.numeric(sub('.', '', fluor$Time))
fluor$Date <- mdy(fluor$Date)
fluor$Time <- fluor$Time * 48 /24

2.2 Visulaize growth:

LD_OTB_all_growth <- ggplot(
 fluor %>%
    filter(
      Treatment %in% c("-Cobalamin", "Replete", "-Thiamine", "-Biotin"),
      Round %in% c(1, 2,3,4)),
  aes(x = Time, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point(size = 3)  + geom_smooth(se=FALSE)+
  # Manual color scale
  scale_color_manual(
    values = c(
      "-Cobalamin" = "#102E50",
      "Replete"    = "#F5C45E",
      "-Thiamine"  = "deeppink4",
      "-Biotin"    = "#667028"
    ),
    labels = c(
      "-Cobalamin" = expression(" —B"[12]),
      "-Thiamine"  = expression("—thiamine"),
      "-Biotin"    = expression("—biotin"),
      "Replete"    = "replete"),
    name = "Treatment") +
 labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  facet_grid(~Round) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.title = element_text(size = 20),
    axis.text  = element_text(size = 20, color = "black"),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 18, color = "black"),
    legend.text  = element_text(size = 18, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

LD_OTB_all_growth

2.3 Calculate maximum RFU values for each round and medium treatment

max_bio <- fluor %>%
  filter(Round %in% c(1, 2, 3, 4)) %>%
  filter(Treatment %in% c("-Cobalamin", "Replete", "-Biotin", "-Thiamine")) %>%
  group_by(Treatment, Round, Rep) %>%
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE), .groups = "drop")

print(max_bio)
## # A tibble: 48 × 4
##    Treatment Round   Rep Max_Chl
##    <chr>     <dbl> <dbl>   <dbl>
##  1 -Biotin       1     1   2911.
##  2 -Biotin       1     2   3220.
##  3 -Biotin       1     3   6239.
##  4 -Biotin       2     1   6771.
##  5 -Biotin       2     2   4230.
##  6 -Biotin       2     3   3492.
##  7 -Biotin       3     1   4288.
##  8 -Biotin       3     2   3237.
##  9 -Biotin       3     3   1767.
## 10 -Biotin       4     1   2049.
## # ℹ 38 more rows

2.3.1 Run an ANOVA on maximum RFU values and Tukey HSD post-hoc test for treatment

aov_max_bio <-  aov(Max_Chl ~ Treatment, data = max_bio)
TukeyHSD(aov_max_bio)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Treatment, data = max_bio)
## 
## $Treatment
##                           diff       lwr       upr     p adj
## -Cobalamin--Biotin   -3365.871 -4602.055 -2129.687 0.0000000
## -Thiamine--Biotin     1961.220   725.036  3197.404 0.0006407
## Replete--Biotin       -757.515 -1993.699   478.669 0.3694455
## -Thiamine--Cobalamin  5327.091  4090.907  6563.275 0.0000000
## Replete--Cobalamin    2608.356  1372.172  3844.540 0.0000068
## Replete--Thiamine    -2718.735 -3954.919 -1482.551 0.0000030
summary(aov_max_bio)
##             Df    Sum Sq  Mean Sq F value   Pr(>F)    
## Treatment    3 174966709 58322236   45.35 1.63e-13 ***
## Residuals   44  56590796  1286154                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.3.2 Prepare data for plotting maximum RFU values:

# Calculate mean and SD
summary_stats <- max_bio %>%
  group_by(Treatment, Round) %>%
  summarise(
    mean_chl = mean(Max_Chl, na.rm = TRUE),
    sd_chl   = sd(Max_Chl, na.rm = TRUE),
    .groups = "drop"
  )

# Make Treatment a factor in both datasets
max_bio$Treatment <- factor(max_bio$Treatment, 
                            levels = c( "-Biotin", "-Cobalamin","-Thiamine", "Replete"))

summary_stats$Treatment <- factor(summary_stats$Treatment, 
                                  levels = c("-Cobalamin", "-Thiamine", "-Biotin", "Replete"))

max_bio$Treatment <- factor(
  max_bio$Treatment,
  levels = c("-Cobalamin", "Replete", "-Thiamine", "-Biotin")
)

summary_stats$Treatment <- factor(
  summary_stats$Treatment,
  levels = c("-Cobalamin", "Replete", "-Thiamine", "-Biotin"))

2.3.3 Plot maximum RFU values:

max_bio_LD <- ggplot(max_bio, aes(x = Treatment, y = Max_Chl, color = Treatment)) +
  
  # Raw points (per round)
  geom_point(size = 3, shape = 16, position = position_jitter(width = 0.1)) +
  
  # Error bars: mean ± SD across replicates
  geom_errorbar(
    data = summary_stats,
    aes(
      x = Treatment,
      y = mean_chl,
      ymin = mean_chl - sd_chl,
      ymax = mean_chl + sd_chl,
      color = Treatment
    ),
    width = 0.3,
    position = position_dodge(width = 0.6),
    inherit.aes = FALSE
  ) +
  
  # Open circles for replicate-averaged means
  geom_point(
    data = summary_stats,
    aes(x = Treatment, y = mean_chl, color = Treatment),
    shape = 21,
    size = 5,
    stroke = 1.2,
    position = position_dodge(width = 0.6),
    inherit.aes = FALSE
  ) +
  
  # Labels
  labs(
    x = "Media Treatment",
    y = expression(paste("Maximum RFU"))
  ) +
  
  # Manual colors
  scale_color_manual(
    values = c(
      "-Cobalamin" = "#102E50",
      "-Thiamine"  = "deeppink4",
      "-Biotin"    = "#667028",
      "Replete"    = "#F5C45E"
    ),
    labels = c(
      "-Cobalamin" = expression("—B"[12]),
      "-Thiamine"  = expression("—Thiamine"),
      "-Biotin"    = expression("—Biotin"),
      "Replete"    = "Replete"
    ),
    name = "Treatment",
    guide = guide_legend(
      override.aes = list(
        shape = 21,
        fill  = "white",
        stroke = 1.2,
        size = 2
      )
    )
  ) +
  
  # X-axis label formatting
  scale_x_discrete(
    labels = c(
      "-Biotin"    = expression("\u2013Biotin"),
      "-Cobalamin" = expression("\u2013B"[12]),
      "-Thiamine"  = expression("\u2013Thiamine"),
      "Replete"    = "Replete"
    )
  ) +
  
  scale_y_continuous(labels = comma) +
  
  # Theme
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 16, color = "black"),
  
    axis.title.y = element_text(size = 22),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 18),
    legend.title = element_text( size = 20, color = "black"),
    legend.text = element_text( size = 16, color = "black"),
    legend.position = "none",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +
  guides(color = "none", fill = "none", shape = "none")+
coord_cartesian(ylim = c(0, 10000)) 

# Print plot
max_bio_LD

2.4 Specific growth rates using all_easylinear

Specific Growth Rates Biotin trial 1 using all_easylinear:

bio1 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(1) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio1, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin1, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(bio_lin1)
## $`1:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -0.1426 -0.3358  0.3733  0.3348  0.1439 -0.2512 -0.1225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.66156    0.21185  17.284 1.19e-05 ***
## x            0.20749    0.02938   7.063  0.00088 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3109 on 5 degrees of freedom
## Multiple R-squared:  0.9089, Adjusted R-squared:  0.8907 
## F-statistic: 49.88 on 1 and 5 DF,  p-value: 0.0008799
## 
## 
## $`1:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -0.2611  0.1266  0.1952  0.1704 -0.1053 -0.1483  0.0224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.77140    0.13104   28.78 9.49e-07 ***
## x            0.21617    0.01817   11.90 7.40e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1923 on 5 degrees of freedom
## Multiple R-squared:  0.9659, Adjusted R-squared:  0.959 
## F-statistic: 141.5 on 1 and 5 DF,  p-value: 7.396e-05
## 
## 
## $`1:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
## -0.018655 -0.066394  0.061972  0.130916 -0.079399 -0.037937  0.009497 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.479439   0.142763   17.37 1.16e-05 ***
## x           0.257843   0.007742   33.30 4.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08194 on 5 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9946 
## F-statistic:  1109 on 1 and 5 DF,  p-value: 4.589e-07
# Extracts only the estimated regression coefficients
coef(bio_lin1)
##                y0    y0_lm     mumax        lag
## 1:-Biotin:1 33.75 38.92188 0.2074863 -0.6871588
## 1:-Biotin:2 33.46 43.44102 0.2161749 -1.2076036
## 1:-Biotin:3 49.69 11.93456 0.2578427  5.5319204
# Returns the model’s R² (proportion of variance explained)
rsquared(bio_lin1)
## 1:-Biotin:1.r2 1:-Biotin:2.r2 1:-Biotin:3.r2 
##      0.9088927      0.9658730      0.9955119

Specific Growth Rates Biotin trial 2 using all_easylinear

bio2 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(2) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio2, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin2, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(bio_lin2)
## $`2:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -1.29805  0.90387  0.62634  0.53091 -0.40115 -0.02684 -0.33508 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.63417    1.14283   1.430  0.21213   
## x            0.37606    0.07849   4.791  0.00492 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8307 on 5 degrees of freedom
## Multiple R-squared:  0.8211, Adjusted R-squared:  0.7854 
## F-statistic: 22.96 on 1 and 5 DF,  p-value: 0.004922
## 
## 
## $`2:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.06137 -0.08479 -0.10821  0.37956 -0.09211  0.27054 -0.30362 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.36851    0.26561  12.682 5.42e-05 ***
## x            0.24466    0.02466   9.921 0.000178 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.261 on 5 degrees of freedom
## Multiple R-squared:  0.9517, Adjusted R-squared:  0.942 
## F-statistic: 98.42 on 1 and 5 DF,  p-value: 0.0001776
## 
## 
## $`2:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.24068 -0.16024 -0.25709 -0.01833  0.16293  0.11462 -0.08257 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.10851    0.16756   24.52 2.10e-06 ***
## x            0.29515    0.01873   15.76 1.87e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1983 on 5 degrees of freedom
## Multiple R-squared:  0.9803, Adjusted R-squared:  0.9763 
## F-statistic: 248.2 on 1 and 5 DF,  p-value: 1.873e-05
# Extracts only the estimated regression coefficients
coef(bio_lin2)
##                 y0     y0_lm     mumax      lag
## 2:-Biotin:1  98.46  5.125203 0.3760631 7.859002
## 2:-Biotin:2  96.62 29.035114 0.2446627 4.914031
## 2:-Biotin:3 107.23 60.856127 0.2951547 1.919209
# Returns the model’s R² (proportion of variance explained)
rsquared(bio_lin2)
## 2:-Biotin:1.r2 2:-Biotin:2.r2 2:-Biotin:3.r2 
##      0.8211458      0.9516543      0.9802552

Specific Growth Rates Biotin trial 3 using all_easylinear:

bio3 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(3) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio3, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin3, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(bio_lin3)
## $`3:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.07856 -0.22592  0.15259  0.26659  0.21258 -0.23432 -0.09296 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.81040    0.35557   7.904 0.000522 ***
## x            0.23963    0.02156  11.115 0.000103 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2282 on 5 degrees of freedom
## Multiple R-squared:  0.9611, Adjusted R-squared:  0.9533 
## F-statistic: 123.5 on 1 and 5 DF,  p-value: 0.0001028
## 
## 
## $`3:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##          1          2          3          4          5          6          7 
##  0.2353339 -0.1112728 -0.0556536 -0.2798767  0.0004855  0.2056337  0.0053499 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.96414    0.23391   12.67 5.44e-05 ***
## x            0.26925    0.01849   14.56 2.76e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1957 on 5 degrees of freedom
## Multiple R-squared:  0.977,  Adjusted R-squared:  0.9723 
## F-statistic:   212 on 1 and 5 DF,  p-value: 2.76e-05
## 
## 
## $`3:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.01417 -0.22670  0.78580  0.49481 -1.90075  0.22234  0.61033 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.56645    2.28724   0.248    0.814
## x            0.18591    0.09401   1.978    0.105
## 
## Residual standard error: 0.9949 on 5 degrees of freedom
## Multiple R-squared:  0.4389, Adjusted R-squared:  0.3267 
## F-statistic: 3.911 on 1 and 5 DF,  p-value: 0.1049
# Extracts only the estimated regression coefficients
coef(bio_lin3)
##                y0     y0_lm     mumax       lag
## 3:-Biotin:1 98.82 16.616595 0.2396343  7.440080
## 3:-Biotin:2 63.40 19.378045 0.2692450  4.402396
## 3:-Biotin:3 82.43  1.761997 0.1859051 20.685290
# Returns the model’s R² (proportion of variance explained)
rsquared(bio_lin3)
## 3:-Biotin:1.r2 3:-Biotin:2.r2 3:-Biotin:3.r2 
##      0.9611016      0.9769567      0.4388915

Specific Growth Rates Biotin trial 4 using all_easylinear

bio4 <- fluor %>%
  filter(Treatment == "-Biotin", Round %in% c(4) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

bio_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, bio4, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(bio_lin4, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(bio_lin4)
## $`4:-Biotin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.30806  0.04318 -0.14397 -0.50617 -0.10602  0.24216  0.16276 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.26650    0.53423   2.371 0.063898 .  
## x            0.23271    0.02897   8.032 0.000484 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3066 on 5 degrees of freedom
## Multiple R-squared:  0.9281, Adjusted R-squared:  0.9137 
## F-statistic: 64.51 on 1 and 5 DF,  p-value: 0.0004838
## 
## 
## $`4:-Biotin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.19446  0.54358 -0.06320 -0.20050 -0.37665  0.05642  0.23480 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.80368    0.59222   3.046  0.02857 * 
## x            0.21685    0.03212   6.752  0.00108 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3399 on 5 degrees of freedom
## Multiple R-squared:  0.9012, Adjusted R-squared:  0.8814 
## F-statistic: 45.58 on 1 and 5 DF,  p-value: 0.001082
## 
## 
## $`4:-Biotin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.02364  0.05129  0.03152 -0.14500 -0.02402 -0.04135  0.10391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.488502   0.120401   20.67 4.91e-06 ***
## x           0.233579   0.008269   28.25 1.04e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08751 on 5 degrees of freedom
## Multiple R-squared:  0.9938, Adjusted R-squared:  0.9925 
## F-statistic: 797.9 on 1 and 5 DF,  p-value: 1.041e-06
# Extracts only the estimated regression coefficients
coef(bio_lin4)
##                y0     y0_lm     mumax       lag
## 4:-Biotin:1 57.52  3.548409 0.2327066 11.970583
## 4:-Biotin:2 67.67  6.071954 0.2168469 11.118270
## 4:-Biotin:3 97.55 12.043218 0.2335787  8.955711
# Returns the model’s R² (proportion of variance explained)
rsquared(bio_lin4)
## 4:-Biotin:1.r2 4:-Biotin:2.r2 4:-Biotin:3.r2 
##      0.9280693      0.9011555      0.9937725

Specific Growth Rates Thiamine trial 1 using all_easylinear:

thia1 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(1) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia1, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin1, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(thia_lin1)
## $`1:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.30173 -0.51432 -0.04376  0.21123  0.09790  0.10677 -0.15955 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4499     0.4668   7.391 0.000713 ***
## x             0.2417     0.0283   8.539 0.000363 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2995 on 5 degrees of freedom
## Multiple R-squared:  0.9358, Adjusted R-squared:  0.923 
## F-statistic: 72.91 on 1 and 5 DF,  p-value: 0.0003627
## 
## 
## $`1:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.31890  0.23129  0.13999 -0.02326  0.25846 -0.25014 -0.03743 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.14430    0.25181   12.49 5.84e-05 ***
## x            0.26660    0.02338   11.40 9.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2474 on 5 degrees of freedom
## Multiple R-squared:  0.963,  Adjusted R-squared:  0.9556 
## F-statistic:   130 on 1 and 5 DF,  p-value: 9.079e-05
## 
## 
## $`1:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.25586  0.36063  0.14809 -0.29385  0.02147 -0.02179  0.04130 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.33255    0.25195   13.23 4.41e-05 ***
## x            0.26006    0.02339   11.12 0.000103 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2476 on 5 degrees of freedom
## Multiple R-squared:  0.9611, Adjusted R-squared:  0.9533 
## F-statistic: 123.6 on 1 and 5 DF,  p-value: 0.0001027
# Extracts only the estimated regression coefficients
coef(thia_lin1)
##                   y0    y0_lm     mumax       lag
## 1:-Thiamine:1 173.20 31.49697 0.2416543 7.0536952
## 1:-Thiamine:2  44.15 23.20350 0.2666038 2.4129054
## 1:-Thiamine:3  33.19 28.00966 0.2600576 0.6525443
# Returns the model’s R² (proportion of variance explained)
rsquared(thia_lin1)
## 1:-Thiamine:1.r2 1:-Thiamine:2.r2 1:-Thiamine:3.r2 
##        0.9358222        0.9629722        0.9611167

Specific Growth Rates Thiamine trial 2 using all_easylinear

thia2 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(2) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia2, h = 9, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin2, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(thia_lin2)
## $`2:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4806 -0.2537 -0.1203  0.2472  0.5625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.63297    0.26430   17.53 4.84e-07 ***
## x            0.24127    0.02348   10.27 1.79e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3638 on 7 degrees of freedom
## Multiple R-squared:  0.9378, Adjusted R-squared:  0.9289 
## F-statistic: 105.5 on 1 and 7 DF,  p-value: 1.79e-05
## 
## 
## $`2:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35944 -0.06839  0.00212  0.08203  0.24602 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.80494    0.11101   43.28 9.17e-10 ***
## x            0.22326    0.01166   19.15 2.64e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1806 on 7 degrees of freedom
## Multiple R-squared:  0.9813, Adjusted R-squared:  0.9786 
## F-statistic: 366.7 on 1 and 7 DF,  p-value: 2.636e-07
## 
## 
## $`2:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42503 -0.12665 -0.00432  0.16254  0.40225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.54004    0.16251   27.94 1.93e-08 ***
## x            0.25614    0.01707   15.01 1.40e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2644 on 7 degrees of freedom
## Multiple R-squared:  0.9699, Adjusted R-squared:  0.9656 
## F-statistic: 225.2 on 1 and 7 DF,  p-value: 1.4e-06
# Extracts only the estimated regression coefficients
coef(thia_lin2)
##                   y0     y0_lm     mumax       lag
## 2:-Thiamine:1 156.08 102.81917 0.2412673 1.7300188
## 2:-Thiamine:2 131.83 122.11205 0.2232630 0.3429772
## 2:-Thiamine:3 140.09  93.69448 0.2561405 1.5704110
# Returns the model’s R² (proportion of variance explained)
rsquared(thia_lin2)
## 2:-Thiamine:1.r2 2:-Thiamine:2.r2 2:-Thiamine:3.r2 
##        0.9378050        0.9812701        0.9698585

Specific Growth Rates Thiamine trial 3 using all_easylinear:

thia3 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(3) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia3, h = 9, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin3, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(thia_lin3)
## $`3:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35661 -0.21370 -0.07224  0.18665  0.48952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.79058    0.36271   4.937  0.00168 ** 
## x            0.25375    0.01937  13.101 3.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3001 on 7 degrees of freedom
## Multiple R-squared:  0.9608, Adjusted R-squared:  0.9552 
## F-statistic: 171.6 on 1 and 7 DF,  p-value: 3.521e-06
## 
## 
## $`3:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49865  0.01803  0.07384  0.10434  0.40557 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.81256    0.35153   5.156  0.00132 ** 
## x            0.31014    0.02091  14.833 1.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3239 on 7 degrees of freedom
## Multiple R-squared:  0.9692, Adjusted R-squared:  0.9648 
## F-statistic:   220 on 1 and 7 DF,  p-value: 1.517e-06
## 
## 
## $`3:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4806 -0.1946 -0.0558  0.3160  0.5117 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.99093    0.40297   4.941  0.00167 ** 
## x            0.25389    0.02397  10.593 1.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3713 on 7 degrees of freedom
## Multiple R-squared:  0.9413, Adjusted R-squared:  0.9329 
## F-statistic: 112.2 on 1 and 7 DF,  p-value: 1.462e-05
# Extracts only the estimated regression coefficients
coef(thia_lin3)
##                  y0    y0_lm     mumax      lag
## 3:-Thiamine:1 49.92 5.992945 0.2537505 8.354027
## 3:-Thiamine:2 72.28 6.126119 0.3101377 7.957711
## 3:-Thiamine:3 57.00 7.322315 0.2538908 8.082708
# Returns the model’s R² (proportion of variance explained)
rsquared(thia_lin3)
## 3:-Thiamine:1.r2 3:-Thiamine:2.r2 3:-Thiamine:3.r2 
##        0.9608132        0.9691648        0.9412777

Specific Growth Rates Thiamine trial 4 using all_easylinear

thia4 <- fluor %>%
  filter(Treatment == "-Thiamine", Round %in% c(4) )%>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

thia_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, thia4, h = 9, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(thia_lin4, log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(thia_lin4)
## $`4:-Thiamine:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34954 -0.14112  0.08196  0.17775  0.22950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.4689     0.1816   19.11 2.68e-07 ***
## x             0.2533     0.0139   18.22 3.71e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2153 on 7 degrees of freedom
## Multiple R-squared:  0.9794, Adjusted R-squared:  0.9764 
## F-statistic:   332 on 1 and 7 DF,  p-value: 3.71e-07
## 
## 
## $`4:-Thiamine:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.269152 -0.100510  0.007589  0.106550  0.281775 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.87449    0.21124   13.61 2.72e-06 ***
## x            0.24309    0.01256   19.35 2.46e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1946 on 7 degrees of freedom
## Multiple R-squared:  0.9816, Adjusted R-squared:  0.979 
## F-statistic: 374.3 on 1 and 7 DF,  p-value: 2.457e-07
## 
## 
## $`4:-Thiamine:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7143  0.0431  0.1993  0.3349  0.4359 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  2.65854    0.76560   3.472   0.0104 * 
## x            0.23859    0.04554   5.239   0.0012 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7055 on 7 degrees of freedom
## Multiple R-squared:  0.7968, Adjusted R-squared:  0.7678 
## F-statistic: 27.45 on 1 and 7 DF,  p-value: 0.0012
# Extracts only the estimated regression coefficients
coef(thia_lin4)
##                   y0    y0_lm     mumax      lag
## 4:-Thiamine:1 161.92 32.10143 0.2532588 6.389518
## 4:-Thiamine:2  76.34 17.71644 0.2430897 6.008910
## 4:-Thiamine:3 109.65 14.27538 0.2385910 8.544988
# Returns the model’s R² (proportion of variance explained)
rsquared(thia_lin4)
## 4:-Thiamine:1.r2 4:-Thiamine:2.r2 4:-Thiamine:3.r2 
##        0.9793538        0.9816432        0.7968191

Specific Growth Rates Replete trial 1 using all_easylinear:

rep1 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(1)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin1 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep1, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin1, log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(rep_lin1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.23658  0.13397  0.19451 -0.01732 -0.04291  0.10941 -0.14107 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.29023    0.23491   14.01 3.34e-05 ***
## x            0.21573    0.01613   13.37 4.19e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1707 on 5 degrees of freedom
## Multiple R-squared:  0.9728, Adjusted R-squared:  0.9674 
## F-statistic: 178.8 on 1 and 5 DF,  p-value: 4.186e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.01471  0.09138 -0.04521 -0.16653  0.07629  0.15920 -0.10042 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.03875    0.17358   17.51 1.12e-05 ***
## x            0.22585    0.01192   18.95 7.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1262 on 5 degrees of freedom
## Multiple R-squared:  0.9863, Adjusted R-squared:  0.9835 
## F-statistic: 358.9 on 1 and 5 DF,  p-value: 7.551e-06
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
##  0.025927  0.009719  0.219121 -0.438641  0.040203  0.154873 -0.011203 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.29711    0.31788   10.37 0.000143 ***
## x            0.22019    0.02183   10.09 0.000164 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.231 on 5 degrees of freedom
## Multiple R-squared:  0.9531, Adjusted R-squared:  0.9438 
## F-statistic: 101.7 on 1 and 5 DF,  p-value: 0.0001641
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(rep_lin1)
##                y0    y0_lm     mumax      lag
## 1:Replete:1 66.01 26.84907 0.2157266 4.169977
## 1:Replete:2 87.82 20.87905 0.2258538 6.360502
## 1:Replete:3 94.03 27.03446 0.2201909 5.661004
# Returns the model R² value (proportion of variance explained)
rsquared(rep_lin1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9727947      0.9862601      0.9531494

Specific Growth Rates Replete trial 2 using all_easylinear:

rep2 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(2)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin2 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep2, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin2, log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(rep_lin2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.18529 -0.03780 -0.16478 -0.18114  0.17936 -0.07894  0.09801 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.37352    0.11553   37.86 2.42e-07 ***
## x            0.18726    0.01602   11.69 8.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1695 on 5 degrees of freedom
## Multiple R-squared:  0.9647, Adjusted R-squared:  0.9576 
## F-statistic: 136.6 on 1 and 5 DF,  p-value: 8.053e-05
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.04068  0.56895 -0.48634 -0.31279 -0.05303  0.38914 -0.06524 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.36733    0.27564  15.844 1.82e-05 ***
## x            0.20689    0.03822   5.412  0.00291 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4045 on 5 degrees of freedom
## Multiple R-squared:  0.8542, Adjusted R-squared:  0.825 
## F-statistic: 29.29 on 1 and 5 DF,  p-value: 0.002913
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.17429  0.11485  0.12945  0.04594 -0.14406  0.10395 -0.07585 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.38836    0.22011   15.39  2.1e-05 ***
## x            0.24033    0.01335   18.01  9.7e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1412 on 5 degrees of freedom
## Multiple R-squared:  0.9848, Adjusted R-squared:  0.9818 
## F-statistic: 324.3 on 1 and 5 DF,  p-value: 9.701e-06
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(rep_lin2)
##                y0    y0_lm     mumax        lag
## 2:Replete:1 95.47 79.32272 0.1872621  0.9894552
## 2:Replete:2 75.69 78.83269 0.2068881 -0.1966365
## 2:Replete:3 51.77 29.61731 0.2403281  2.3237062
# Returns the model R² value (proportion of variance explained)
rsquared(rep_lin2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9646965      0.8542034      0.9848147

Specific Growth Rates Replete trial 3 using all_easylinear

rep3 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(3)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin3 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep3, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin3, log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(rep_lin3)
## $`3:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5         6         7 
##  0.056275  0.286364 -0.388386  0.001474 -0.003777 -0.212798  0.260848 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.26798    0.31594   7.179 0.000816 ***
## x            0.26516    0.02498  10.616 0.000128 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2643 on 5 degrees of freedom
## Multiple R-squared:  0.9575, Adjusted R-squared:  0.949 
## F-statistic: 112.7 on 1 and 5 DF,  p-value: 0.0001282
## 
## 
## $`3:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.01109 -0.07478 -0.04567  0.22757 -0.09197 -0.00871 -0.01753 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.69027    0.20383   8.292 0.000416 ***
## x            0.25682    0.01105  23.232 2.75e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.117 on 5 degrees of freedom
## Multiple R-squared:  0.9908, Adjusted R-squared:  0.989 
## F-statistic: 539.7 on 1 and 5 DF,  p-value: 2.75e-06
## 
## 
## $`3:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.37438 -0.38671 -0.12656 -0.27670  0.47488  0.07385 -0.13314 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.09308    0.42459    4.93 0.004361 ** 
## x            0.30144    0.03357    8.98 0.000286 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3552 on 5 degrees of freedom
## Multiple R-squared:  0.9416, Adjusted R-squared:  0.9299 
## F-statistic: 80.65 on 1 and 5 DF,  p-value: 0.0002856
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(rep_lin3)
##                y0    y0_lm     mumax       lag
## 3:Replete:1 46.52 9.659839 0.2651609  5.928118
## 3:Replete:2 77.27 5.420959 0.2568156 10.346075
## 3:Replete:3 61.72 8.109824 0.3014418  6.732749
# Returns the model R² value (proportion of variance explained)
rsquared(rep_lin3)
## 3:Replete:1.r2 3:Replete:2.r2 3:Replete:3.r2 
##      0.9575206      0.9908211      0.9416198

Specific Growth Rates Replete trial 4 using all_easylinear:

rep4 <- fluor %>%
  filter(Treatment == "Replete", Round %in% c(4)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

rep_lin4 <- all_easylinear(Chl ~ Time|Round+Treatment+Rep, rep4, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(rep_lin4, log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(rep_lin4)
## $`4:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.09195  0.19036  0.07953 -0.20968 -0.15320  0.21726 -0.03232 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.00947    0.25121   11.98 7.15e-05 ***
## x            0.21058    0.01725   12.21 6.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1826 on 5 degrees of freedom
## Multiple R-squared:  0.9675, Adjusted R-squared:  0.961 
## F-statistic:   149 on 1 and 5 DF,  p-value: 6.529e-05
## 
## 
## $`4:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.11124  0.04917 -0.34079  0.14249 -0.10230  0.22701 -0.08682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.08172    0.28910   10.66 0.000126 ***
## x            0.20366    0.01986   10.26 0.000151 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2101 on 5 degrees of freedom
## Multiple R-squared:  0.9546, Adjusted R-squared:  0.9456 
## F-statistic: 105.2 on 1 and 5 DF,  p-value: 0.0001513
## 
## 
## $`4:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.14336  0.08611 -0.35495  0.21552 -0.30373  0.09004  0.12366 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.07289    0.43805   4.732 0.005187 ** 
## x            0.17495    0.02376   7.364 0.000725 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2514 on 5 degrees of freedom
## Multiple R-squared:  0.9156, Adjusted R-squared:  0.8987 
## F-statistic: 54.23 on 1 and 5 DF,  p-value: 0.0007253
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(rep_lin4)
##                 y0     y0_lm     mumax       lag
## 4:Replete:1  78.16 20.276716 0.2105809  6.407442
## 4:Replete:2 157.37 21.795906 0.2036571  9.706892
## 4:Replete:3  87.34  7.947773 0.1749477 13.700762
# Returns the model R² value (proportion of variance explained)
rsquared(rep_lin4)
## 4:Replete:1.r2 4:Replete:2.r2 4:Replete:3.r2 
##      0.9675256      0.9546300      0.9155837

2.4.1 Combine all specific growth rate values from each treatment:

extract_results <- function(model, treatment_name) {
  coefs <- coef(model)
  r2vals <- rsquared(model)
  
  df <- as.data.frame(coefs) %>%
    tibble::rownames_to_column("Group") %>%
    mutate(
      Treatment = treatment_name,
      r2 = r2vals
    ) %>%
    # Split "Group" into Round and Rep (Group looks like "1:-Thiamine:2")
    tidyr::separate(Group, into = c("Round", "Treatment_label", "Rep"), sep = ":", remove = FALSE) %>%
    mutate(Round = as.integer(Round))
  
  return(df)
}

# --- Extract Biotin rounds ---
df_bio1 <- extract_results(bio_lin1, "-Biotin")
df_bio2 <- extract_results(bio_lin2, "-Biotin")
df_bio3 <- extract_results(bio_lin3, "-Biotin")
df_bio4 <- extract_results(bio_lin4, "-Biotin")
df_bio <- bind_rows(df_bio1, df_bio2, df_bio3, df_bio4)

# --- Extract Thiamine rounds ---
df_thia1 <- extract_results(thia_lin1, "-Thiamine")
df_thia2 <- extract_results(thia_lin2, "-Thiamine")
df_thia3 <- extract_results(thia_lin3, "-Thiamine")
df_thia4 <- extract_results(thia_lin4, "-Thiamine")
df_thia <- bind_rows(df_thia1, df_thia2, df_thia3, df_thia4)

# --- Extract Replete rounds ---
df_rep1 <- extract_results(rep_lin1, "Replete")
df_rep2 <- extract_results(rep_lin2, "Replete")
df_rep3 <- extract_results(rep_lin3, "Replete")
df_rep4 <- extract_results(rep_lin4, "Replete")
df_rep <- bind_rows(df_rep1, df_rep2, df_rep3, df_rep4)

# --- Combine all treatments ---
all_lin_results <- bind_rows(df_bio, df_thia, df_rep)

# Optional: ensure Round is treated as a factor
all_lin_results$Round <- factor(all_lin_results$Round)

2.4.2 Run an ANOVA on specific growth rate values and Tukey HSD post-hoc test for treatment:

all_lin_results$Round <- factor(all_lin_results$Round)
linear <- aov(mumax ~ Treatment , data = all_lin_results)
summary(linear)
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Treatment    2 0.00517 0.002586   1.861  0.172
## Residuals   33 0.04587 0.001390
TukeyHSD(linear)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mumax ~ Treatment, data = all_lin_results)
## 
## $Treatment
##                           diff         lwr         upr     p adj
## -Thiamine--Biotin  0.005533672 -0.03181359 0.042880937 0.9298833
## Replete--Biotin   -0.022203946 -0.05955121 0.015143319 0.3234309
## Replete--Thiamine -0.027737618 -0.06508488 0.009609647 0.1780243

2.4.3 Plot specific growth rates for Pyro grown in -Thiamine, -Biotin, and Replete:

all_lin_results$Treatment <- factor(
  all_lin_results$Treatment,
  levels = c("-Cobalamin",  "-Thiamine", "-Biotin","Replete")
)

mu <- ggplot(all_lin_results, aes(x = Treatment, y = mumax, color = Treatment)) +
  # Raw points
  geom_point(size = 3, shape = 16, position = position_jitter(width = 0.1)) +
  # Mean ± SE
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  scale_color_manual(values = c(
    "-Cobalamin" = "#102E50",
    "Replete"    = "#F5C45E",
    "-Thiamine"  = "deeppink4",
    "-Biotin"    = "#667028"
  )) +
  scale_x_discrete(labels = c(
    "-Cobalamin" = expression("\u2212B"[12]),
    "Replete"    = "Replete",
    "-Thiamine"  = expression("\u2212Thiamine"),
    "-Biotin"    = expression("\u2212Biotin")
  )) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text( size = 20),
    axis.text = element_text( size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text( size = 20, color = "black"),
    legend.title = element_blank(),
    legend.text = element_blank(),
    legend.key = element_blank()
  ) +
    labs(
      x = "Media Treatment",
      y = expression(paste("(", mu, ") RFU" * "\u00B7" * "day"^-1, " "))
    ) + guides(color = "none", fill = "none", shape = "none")

mu

2.5 Figure 3

bo <- 
  LD_OTB_all_growth / 
  (max_bio_LD +
   mu ) +
  plot_layout(guides = "collect") &
  theme(
    legend.position = "bottom",
    panel.spacing = unit(.1, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    axis.text.x = element_text(size = 16),
    
  )
bo

Fig. 3 Growth responses of P. bahamense when grown in replete L1/4 media and in media missing vitamin B12, thiamine, or biotin. A In vivo Chlorophyll a fluorescence (Chl-a) was used to track P. bahamense growth across four experimental rounds. The replicates growing in replete media and in media missing thiamine or biotin were serially transferred between each round. Since replicates in media without B12 failed to grow, the replete replicates were used to restart the cultures in media missing B12 in each round. B Maximum Chl-a values in each replicate by media formulation. All data are plotted with color representing the media treatment. Means are plotted as black open circles, and error bars represent one standard deviation of the mean. ANOVA results indicate media had a significant effect on maximum Chl-a (p = 0.0000, F = 45.35). Post-hoc multiple comparisons using the Tukey HSD test indicate that the maximum biomasses achieved in replicates grown in media missing B12 were significantly lower than in replicates grown in replete media (p-adj = 0.0001), media missing thiamine (p-adj = 0.0001), and media missing biotin (p-adj = 0.0001). Replicates grown in media missing thiamine reached significantly higher Chl-a values than replicates grown in replete media (p-adj = 0.0001) or in media without biotin (p-adj = 0.001). The maximum Chl-a values for replicates grown in media missing biotin were not significantly different from replicates in replete media (p-adj = 0.369). Results from comparisons to replete are shown on the plot with the following significance codes: p-adj > 0.05; ‘n.s.’ and p-adj < 0.001; ’*’. C** Specific growth rate (µ) of each replicate in –thiamine, –biotin, and replete L1/4 media. Growth rates could not be calculated for replicates in L1/4 –B12 due to low and variable Chl-a values. Means are plotted as black open circles and error bars represent one standard deviation of the mean. Excluding the –B12 treatment, media treatment did not significantly affect specific growth rate (one-way ANOVA; p = 0.172, F = 1.861).

3 Evaluating the B12 requirements of P. bahamense

Chlorophyll-a fluorescence (RFU) data from P. bahamense cultures were used to visualize growth patterns, quantify maximum biomass, and estimate specific growth rates under media conditions deficient in media lacking one targeted B vitamin. Specific growth rates were computed using the “alleasy_linear” function of the “growthrates” package in R Studio (Hall et al., 2014; RStudio Team, 2020).

3.1 Import and Parse the Dataset:

b12 <-as.data.frame(read_csv("b12_levels.csv")) 
# Remove periods and convert to numeric
b12$Time <- as.numeric(sub("T", "", b12$Time))
# Convert Date column to Date type
b12$Date <- mdy(b12$Date)
b12$Treatment <- unlist(b12$Treatment)
#covert timepoints to hours
b12$Time <- b12$Time * 48    # hours
b12$Time <- b12$Time / 24

# put treatment level in desired order
b12 <- b12 %>%
  mutate(Treatment = factor(Treatment, levels = c(
    "369pM",    # Place the levels in desired order
    "36pM",
    "3.69pM",
    "0.369pM",
    "0.0369pM"
  )))

3.2 Prepare data for plotting:

levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")
b12$Treatment <- factor(b12$Treatment, levels = levels_order)

# Named color palette
b12_palette <- c(
  "0.0369pM" = "lightblue",
  "0.369pM"  = "#99c2e6",
  "3.69pM"   = "#66a3d2",
  "36pM"     = "#3380bf",
  "369pM"    = "#004080"
)

3.3 Visulaize growth

growth <- ggplot(b12, aes(x = Time, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point() +
  geom_smooth(aes(color = Treatment), se = FALSE) +
  facet_grid(Round ~ Treatment) +
  labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  
 
   scale_y_continuous(labels = scales::comma) +


  scale_color_manual(values = b12_palette) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 12, color = "black"),
    axis.title.y = element_text(size = 20, margin = margin(t = 10)),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 16),
    strip.text = element_text( size = 16, color = "black"),
    legend.title = element_text( size = 16, color = "black"),
    legend.text = element_text( size = 16, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt"),
    panel.spacing = unit(.15, "cm")
  )

growth

3.4 Calculate and visualize maximum RFU values for each round and medium treatment:

3.4.1 Calculate maximum RFU:

max_bio <- b12  %>%
  filter(Treatment %in% c(Treatment, levels = c("369pM", "36pM", "3.69pM", "0.369pM", "0.0369pM"))) %>%
  group_by(Treatment, Bio_Rep) %>% 
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE)) %>%
  ungroup()

print(max_bio)
## # A tibble: 15 × 3
##    Treatment Bio_Rep Max_Chl
##    <fct>       <dbl>   <dbl>
##  1 0.0369pM        1    842.
##  2 0.0369pM        2    541.
##  3 0.0369pM        3    928.
##  4 0.369pM         1   1968.
##  5 0.369pM         2   1897.
##  6 0.369pM         3   1866.
##  7 3.69pM          1   1923.
##  8 3.69pM          2   3929.
##  9 3.69pM          3   2146.
## 10 36pM            1   5580.
## 11 36pM            2   6041.
## 12 36pM            3   4801.
## 13 369pM           1   5287.
## 14 369pM           2   4424.
## 15 369pM           3   5539.

3.4.2 Run an ANOVA followed by a Tukey HSD test for the treatment factor:

result<- aov(Max_Chl~Treatment, max_bio)
summary(result)
##             Df   Sum Sq  Mean Sq F value   Pr(>F)    
## Treatment    4 49705730 12426432   31.27 1.25e-05 ***
## Residuals   10  3973244   397324                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(result)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Treatment, data = max_bio)
## 
## $Treatment
##                       diff        lwr      upr     p adj
## 0.369pM-0.0369pM 1139.8033  -554.0107 2833.617 0.2494204
## 3.69pM-0.0369pM  1895.5967   201.7826 3589.411 0.0272120
## 36pM-0.0369pM    4703.7467  3009.9326 6397.561 0.0000277
## 369pM-0.0369pM   4312.9433  2619.1293 6006.757 0.0000597
## 3.69pM-0.369pM    755.7933  -938.0207 2449.607 0.6023509
## 36pM-0.369pM     3563.9433  1870.1293 5257.757 0.0003042
## 369pM-0.369pM    3173.1400  1479.3259 4866.954 0.0007806
## 36pM-3.69pM      2808.1500  1114.3359 4501.964 0.0020036
## 369pM-3.69pM     2417.3467   723.5326 4111.161 0.0058828
## 369pM-36pM       -390.8033 -2084.6174 1303.011 0.9366508

3.4.3 Plot maximum biomass:

# Define factor levels (low to high B12)
levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")
max_bio$Treatment <- factor(trimws(max_bio$Treatment), levels = levels_order)

# Calculate summary statistics
summary_stats <- max_bio %>%
  group_by(Treatment) %>%
  summarise(
    mean_chl = mean(Max_Chl, na.rm = TRUE),
    sd_chl   = sd(Max_Chl, na.rm = TRUE),
    .groups = "drop"
  )
summary_stats$Treatment <- factor(summary_stats$Treatment, levels = levels_order)

# Define color palette
b12_palette <- c(
  "0.0369pM" = "lightblue",
  "0.369pM"  = "#99c2e6",
  "3.69pM"   = "#66a3d2",
  "36pM"     = "#3380bf",
  "369pM"    = "#004080"
)


max_bio_plot <- ggplot(summary_stats, aes(x = Treatment)) +
  
  # Replicate points
  geom_point(
    data = max_bio,
    aes(y = Max_Chl, color = Treatment),
    size = 3,
    alpha = 1,
    position = position_jitter(width = 0.1)
  ) +
  
 stat_summary(
  data = max_bio,
  aes(y = Max_Chl, color = Treatment),
  fun = mean,
  geom = "point",
  shape = 1,
  size = 5,
  stroke = 1.2
) +

stat_summary(
  data = max_bio,
  aes(y = Max_Chl, color = Treatment),
  fun.data = mean_sdl,
  fun.args = list(mult = 1),
  geom = "errorbar",
  width = 0.2
)+
  
  # Keep your palette but remove the scale (so no legend)
  scale_color_manual(values = b12_palette, guide = "none") +
  
  scale_y_continuous(labels = scales::comma) +
  
  coord_cartesian(ylim = c(0, 10000)) +
  
  ylab(expression(paste("Maximum RFU"))) +
  xlab(expression(B[12]~Concentration~"(pM)")) +
  
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 20, color = "black", angle = 45, hjust = 1),
    axis.title.y = element_text(size = 20),
    axis.title.x = element_text(size = 20),
    legend.position = "none",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

max_bio_plot

3.5 Specific growth rates using all_easylinear:

3.5.1 Remove outliers:

# ️⃣ Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(b12, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
b12_merged <- merge(b12, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

#  Keep all T0 points, and remove any post-T0 points where replicate < its T0
b12_clean <- subset(
  b12_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

Specific growth rate with all easy linear for 0.0369pM:

pM_1 <- b12_clean %>%
  filter(
    Treatment %in% c("0.0369pM"),
    Round %in% c(1)
  ) %>%                     
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu1 <- all_easylinear(Chl ~ Time | Round + Treatment + Bio_Rep, pM_1, h =4, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(mu1, log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(mu1)
## $`1:0.0369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.21354 -0.37851  0.11642  0.04856 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  4.72824    0.52579   8.993   0.0121 *
## x            0.19590    0.07155   2.738   0.1115  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.32 on 2 degrees of freedom
## Multiple R-squared:  0.7894, Adjusted R-squared:  0.6841 
## F-statistic: 7.496 on 1 and 2 DF,  p-value: 0.1115
## 
## 
## $`1:0.0369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.10998  0.29391 -0.05595 -0.12799 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.39799    0.22595  23.890  0.00175 **
## x            0.07310    0.02382   3.069  0.09177 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2429 on 2 degrees of freedom
## Multiple R-squared:  0.8249, Adjusted R-squared:  0.7373 
## F-statistic: 9.421 on 1 and 2 DF,  p-value: 0.09177
## 
## 
## $`1:0.0369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.29094 -0.55228  0.23175  0.02959 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   4.5299     0.5772   7.848   0.0159 *
## x             0.2125     0.1054   2.017   0.1813  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4713 on 2 degrees of freedom
## Multiple R-squared:  0.6703, Adjusted R-squared:  0.5055 
## F-statistic: 4.067 on 1 and 2 DF,  p-value: 0.1813
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(mu1)
##                  y0     y0_lm      mumax        lag
## 1:0.0369pM:1 145.02 113.09593 0.19590188  1.2691828
## 1:0.0369pM:2 197.95 220.96199 0.07310236 -1.5044141
## 1:0.0369pM:3  97.44  92.75275 0.21252290  0.2319725
# Returns the model R² value (proportion of variance explained)
rsquared(mu1)
## 1:0.0369pM:1.r2 1:0.0369pM:2.r2 1:0.0369pM:3.r2 
##       0.7893916       0.8248814       0.6703408
#not significant slopes

Specific growth rate with all easy linear for 0.369pM:

pM_2 <- b12_clean %>%
  filter(Treatment %in% c("0.369pM"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu2 <-all_easylinear(Chl ~ Time|Treatment+Bio_Rep, pM_2, h = 4, quota = 1)


par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(mu2, log = "y")   # line only

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(mu2)
## $`0.369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.03871  0.14944 -0.14406  0.03332 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.72896    0.22995  20.566  0.00236 **
## x            0.23253    0.02555   9.101  0.01186 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1512 on 2 degrees of freedom
## Multiple R-squared:  0.9764, Adjusted R-squared:  0.9646 
## F-statistic: 82.83 on 1 and 2 DF,  p-value: 0.01186
## 
## 
## $`0.369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.03225 -0.15070  0.39814 -0.21519 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  5.12741    0.70102   7.314   0.0182 *
## x            0.20224    0.07559   2.675   0.1159  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3381 on 2 degrees of freedom
## Multiple R-squared:  0.7816, Adjusted R-squared:  0.6724 
## F-statistic: 7.157 on 1 and 2 DF,  p-value: 0.1159
## 
## 
## $`0.369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.26163  0.40004 -0.01518 -0.12322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  4.21953    0.72414   5.827   0.0282 *
## x            0.28628    0.07809   3.666   0.0670 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3492 on 2 degrees of freedom
## Multiple R-squared:  0.8705, Adjusted R-squared:  0.8057 
## F-statistic: 13.44 on 1 and 2 DF,  p-value: 0.06701
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(mu2)
##               y0     y0_lm     mumax      lag
## 0.369pM:1 248.73 113.17754 0.2325283 3.386298
## 0.369pM:2 375.27 168.57963 0.2022351 3.956966
## 0.369pM:3 182.06  68.00143 0.2862760 3.440063
# Returns the model R² value (proportion of variance explained)
rsquared(mu2)
## 0.369pM:1.r2 0.369pM:2.r2 0.369pM:3.r2 
##    0.9764234    0.7815961    0.8704739
#not significant slopes

Specific growth rate with all easy linear for 3.69pM:

pM_3 <- b12_clean %>%
  filter(Treatment %in% c("3.69pM"),
         Round %in% c(1),
         Time >= 0) %>%  
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


mu3 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_3, h =5, quota = 1)


par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(mu3 , log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(mu3)
## $`1:3.69pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.03126 -0.21999  0.32699  0.13101 -0.20676 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.60153    0.36017  12.776  0.00103 **
## x            0.24850    0.04245   5.854  0.00994 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2685 on 3 degrees of freedom
## Multiple R-squared:  0.9195, Adjusted R-squared:  0.8927 
## F-statistic: 34.28 on 1 and 3 DF,  p-value: 0.009935
## 
## 
## $`1:3.69pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
##  0.100883 -0.081917 -0.005946 -0.145891  0.132871 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.40262    0.10570   51.11 1.65e-05 ***
## x            0.21586    0.02158   10.01  0.00213 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1365 on 3 degrees of freedom
## Multiple R-squared:  0.9709, Adjusted R-squared:  0.9612 
## F-statistic: 100.1 on 1 and 3 DF,  p-value: 0.002126
## 
## 
## $`1:3.69pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
## -0.076428  0.082299  0.033975 -0.009133 -0.030712 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8058     0.1154   41.65 3.05e-05 ***
## x             0.2069     0.0111   18.63 0.000337 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07023 on 3 degrees of freedom
## Multiple R-squared:  0.9914, Adjusted R-squared:  0.9886 
## F-statistic: 347.1 on 1 and 3 DF,  p-value: 0.0003375
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(mu3)
##                y0     y0_lm     mumax       lag
## 1:3.69pM:1 142.16  99.63649 0.2485048 1.4302532
## 1:3.69pM:2 245.55 221.98664 0.2158566 0.4673628
## 1:3.69pM:3 271.24 122.21605 0.2068770 3.8535635
# Returns the model R² value (proportion of variance explained)
rsquared(mu3)
## 1:3.69pM:1.r2 1:3.69pM:2.r2 1:3.69pM:3.r2 
##     0.9195173     0.9708995     0.9914315

Specific growth rate with all easy linear for 36pM:

pM_4 <- b12_clean %>% 
  filter(Treatment %in% c("36pM"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


mu4 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_4, h = 5, quota = 1)


par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(mu4 , log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(mu4)
## $`1:36pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.02088 -0.04810  0.12932 -0.19786  0.09576 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.69402    0.24700   19.00 0.000318 ***
## x            0.23918    0.02377   10.06 0.002090 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1503 on 3 degrees of freedom
## Multiple R-squared:  0.9712, Adjusted R-squared:  0.9616 
## F-statistic: 101.3 on 1 and 3 DF,  p-value: 0.00209
## 
## 
## $`1:36pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5 
## -0.5827  0.4830  0.2809  0.1390 -0.3201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.92981    0.81583   6.043  0.00909 **
## x            0.18334    0.06122   2.995  0.05791 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5078 on 3 degrees of freedom
## Multiple R-squared:  0.7494, Adjusted R-squared:  0.6658 
## F-statistic: 8.969 on 1 and 3 DF,  p-value: 0.05791
## 
## 
## $`1:36pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.05819  0.08475 -0.32792  0.16882  0.01616 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.22326    0.36356  11.617  0.00137 **
## x            0.27122    0.03498   7.753  0.00446 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2213 on 3 degrees of freedom
## Multiple R-squared:  0.9525, Adjusted R-squared:  0.9366 
## F-statistic: 60.11 on 1 and 3 DF,  p-value: 0.004463
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(mu4)
##              y0     y0_lm     mumax      lag
## 1:36pM:1 330.46 109.29127 0.2391784 4.626125
## 1:36pM:2 205.58 138.35322 0.1833369 2.160096
## 1:36pM:3 104.26  68.25567 0.2712184 1.561941
# Returns the model R² value (proportion of variance explained)
rsquared(mu4)
## 1:36pM:1.r2 1:36pM:2.r2 1:36pM:3.r2 
##   0.9712271   0.7493526   0.9524610

Specific growth rate with all easy linear for 369pM:

pM_5 <- b12_clean %>% 
  filter(Treatment %in% c("369pM"),
         Round %in% c(1),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

mu5 <-all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, pM_5, h = 5, quota = 1)


par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(mu5 , log = "y")

# Shows the full regression output: coefficients, SEs, t-tests, p-values, and model fit statistics
summary(mu5)
## $`1:369pM:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.08633  0.17016 -0.04861  0.01842 -0.05363 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.74188    0.15680   30.24 7.94e-05 ***
## x            0.25385    0.01538   16.51 0.000484 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1183 on 3 degrees of freedom
## Multiple R-squared:  0.9891, Adjusted R-squared:  0.9855 
## F-statistic: 272.6 on 1 and 3 DF,  p-value: 0.0004836
## 
## 
## $`1:369pM:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##         1         2         3         4         5 
##  0.221300 -0.358207  0.001927 -0.035733  0.170713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6394     0.3488  13.303 0.000918 ***
## x             0.2400     0.0342   7.017 0.005946 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2631 on 3 degrees of freedom
## Multiple R-squared:  0.9426, Adjusted R-squared:  0.9234 
## F-statistic: 49.23 on 1 and 3 DF,  p-value: 0.005946
## 
## 
## $`1:369pM:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.26764  0.51741 -0.14472 -0.19223  0.08718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  4.27559    0.49285   8.675  0.00322 **
## x            0.25738    0.05808   4.431  0.02135 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3673 on 3 degrees of freedom
## Multiple R-squared:  0.8675, Adjusted R-squared:  0.8233 
## F-statistic: 19.64 on 1 and 3 DF,  p-value: 0.02135
# Extracts only the estimated regression coefficients (intercept + predictors)
coef(mu5)
##               y0     y0_lm     mumax      lag
## 1:369pM:1 286.00 114.64947 0.2538508 3.600983
## 1:369pM:2 321.18 103.48174 0.2399582 4.720016
## 1:369pM:3 106.96  71.92225 0.2573820 1.541947
# Returns the model R² value (proportion of variance explained)
rsquared(mu5)
## 1:369pM:1.r2 1:369pM:2.r2 1:369pM:3.r2 
##    0.9891142    0.9425652    0.8674686

3.5.2 Make dataframe with fits:

  make_df <- function(model, round_num, experiment_name) {
  coefs <- coef(model)
  
  if (is.matrix(coefs)) {
    df <- as.data.frame(coefs) %>%
      tibble::rownames_to_column("Group") %>%
      mutate(
        Round = round_num,
        Experiment = experiment_name,
        r2 = rsquared(model)
      )
  } else {
    df <- data.frame(
      Group = names(coefs),
      Estimate = as.numeric(coefs),
      Round = round_num,
      Experiment = experiment_name,
      r2 = rep(rsquared(model), length(coefs))
    )
  }
  
  return(df)
}

# Example usage for your single fits
df1 <- make_df(mu1, 1, "0.0369pM")
df2 <- make_df(mu2, 1, "0.369pM")
df3 <- make_df(mu3, 1, "3.69pM")
df4 <- make_df(mu4, 1, "36pM")
df5 <- make_df(mu5, 1, "369pM")

# Combine all results
all_results <- bind_rows(df1, df2, df3, df4, df5)

# Check
all_results

3.5.3 ANOVA and post-hoc multiple comparisons of mumax among varying media treatments:

aov_result <- aov(mumax ~ Experiment, data = all_results)
summary(aov_result)
##             Df  Sum Sq  Mean Sq F value Pr(>F)
## Experiment   4 0.01503 0.003757    1.85  0.196
## Residuals   10 0.02031 0.002031
# Post-hoc Tukey test
TukeyHSD(aov_result)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = mumax ~ Experiment, data = all_results)
## 
## $Experiment
##                          diff         lwr       upr     p adj
## 0.369pM-0.0369pM  0.079837439 -0.04125361 0.2009285 0.2651098
## 3.69pM-0.0369pM   0.063237082 -0.05785396 0.1843281 0.4652490
## 369pM-0.0369pM    0.089887949 -0.03120310 0.2109790 0.1807721
## 36pM-0.0369pM     0.070735492 -0.05035555 0.1918265 0.3656496
## 3.69pM-0.369pM   -0.016600357 -0.13769140 0.1044907 0.9901101
## 369pM-0.369pM     0.010050511 -0.11104054 0.1311416 0.9985578
## 36pM-0.369pM     -0.009101947 -0.13019299 0.1119891 0.9990215
## 369pM-3.69pM      0.026650867 -0.09444018 0.1477419 0.9459297
## 36pM-3.69pM       0.007498410 -0.11359264 0.1285895 0.9995435
## 36pM-369pM       -0.019152457 -0.14024350 0.1019386 0.9831961

3.5.4 T-test between 0.0369pM and all other concentrations:

all_results$Group <- ifelse(all_results$Experiment == "0.0369pM", 
                            "A", 
                            "B")
all_results$Group <- factor(all_results$Group)

t.test(mumax~ Group, data = all_results)
## 
##  Welch Two Sample t-test
## 
## data:  mumax by Group
## t = -1.6943, df = 2.1577, p-value = 0.2231
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
##  -0.2558391  0.1039901
## sample estimates:
## mean in group A mean in group B 
##       0.1605090       0.2364335

3.5.5 Plot mumax with all easy linear data:

# Define the order of treatments from least to greatest
levels_order <- c("0.0369pM", "0.369pM", "3.69pM", "36pM", "369pM")

# Convert Treatment column to a factor with this order
all_results$Experiment <- factor(all_results$Experiment, levels = levels_order)

# calculate standard deviation and mean 
b12_summary <- all_results %>%
  group_by(Experiment) %>%
  summarise(
    mean_mumax = mean(mumax, na.rm = TRUE),
    se_mumax   = sd(mumax, na.rm = TRUE)/sqrt(n()),
    .groups = "drop"
  )

mu_b12 <- ggplot(all_results, aes(x = Experiment, y = mumax, color = Experiment)) +
  
  # Raw replicate points (jittered)
  geom_point(size = 3, shape = 16, 
             position = position_jitter(width = 0.1), alpha = 1) +
  
  # Mean per Experiment+Round (open circle)
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2)  +
  
  # Error bars (standard error)
  stat_summary(fun.data = mean_se, geom = "errorbar",
               width = 0.2) +
  
  # Color palette — WITH legend removed
  scale_color_manual(
    values = c(
      "0.0369pM" = "lightblue",
      "0.369pM"  = "#99c2e6",
      "3.69pM"   = "#66a3d2",
      "36pM"     = "#3380bf",
      "369pM"    = "#004080"
    ),
    guide = "none"   # ← THIS removes the legend safely
  ) +

  # Safe visible range
  coord_cartesian(ylim = c(0, 0.4)) +
  
  # Labels + theme
  labs(
    x = expression(B[12]~Concentration~"(pM)"),
    y = expression(paste(mu, " (RFU" %*% "day"^-1, ")"))
  ) +
  
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text(size = 20),
    axis.text = element_text(size = 20, color = "black"),
    axis.text.x = element_text(size = 20, color = "black", angle = 45, hjust = 1),
    strip.text = element_text(size = 20, color = "black"),

    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  )

mu_b12

3.6 Figure 4

# Ensure BOTH subplots have no legend inside patchwork
max_bio_nl <- max_bio_plot + theme(legend.position = "none")
mu_nl      <- mu_b12 + theme(legend.position = "none")

# Build the layout
bo <- 
  (
    growth /                
    (max_bio_nl | mu_nl)    
  ) +
  plot_layout(
    widths = c(4, 2),        # your original width ratio
    guides = "collect"        # ← THIS is what actually gathers one legend
  ) &
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    panel.spacing = unit(0.4, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    legend.ticks = element_blank(),
    axis.text.x = element_text(size = 16)
  )

bo

Fig. 4 Growth responses of P. bahamense when grown in L1/4 media with a 10% dilution series of B12 concentrations. A In vivo Chlorophyll a fluorescence (Chl-a) was used to track P. bahamense growth in L1/4 media with five concentrations of B12: 0.0369 pM, 0.369 pM, 36 pM, and 369 pM. B Maximum Chl-a values in each replicate by media formulation. Means are plotted as black open circles, and error bars represent one standard deviation of the mean. One-way ANOVA results indicate media had a significant effect on maximum Chl-a (p = 0.0000, F = 31.27 ). Post-hoc multiple comparisons using the Tukey HSD test indicate that the maximum Chl-a achieved in replicates grown in media with the lowest concentration of B12 (0.0369 pM) was significantly lower than in media with 3.69, 36, and 369 pM B12 (p-adj = 0.025, 0.0000, 0.0001, respectively). There were no significant differences in maximum Chl-a achieved between successive B12 concentrations except between 3.69 and 36 pM (p-adj = 0.002). Results from key comparisons are shown on the plot with the following significance codes: p-adj > 0.05; ‘n.s.’, p-adj < 0.05; ’’, and p-adj < 0.001; ’’. C** Specific growth rates (µ; RFU day–1) for each replicate by media formulation. All data are plotted with color representing the media treatment. Means are plotted as black open circles, and error bars represent one standard deviation of the mean. ANOVA results indicate media did not have a significant effect on growth rate (p = 0.196, F = 1.85). However, the growth rates for the replicates grown in media with 0.0369 pM B12 were poor fits (p = 0.1813, R2 < 0.85). The poor fits could have led to overestimated growth rate, obscuring the depressed growth rate in 0.0369 pM B12 and skewing the results of statistical tests.

4 Evaluation of bacterial community potential to supplement B12 and support P. bahamense blooms

Low-diveristy P. bahamense culture was co-cultured with bacteria recruited from OTB (Old Tampa Bay) water, Indian River Lagoon (IRL) culture, and IRL surface seawater and grown in replete media and media missing B12

4.1 Import data:

OTB water:

OTB <- as.data.frame(read_csv("OTB_1006.csv"))

OTB$Time <- as.numeric(sub('.', '', OTB$Time))
OTB$Date <- ymd(OTB$Date)
OTB$Time <- OTB$Time * 48

OTB %<>%
  filter(Treatment %in% c("-Cobalamin", "Replete")) %>%
  mutate(Treatment = factor(Treatment, levels = c("-Cobalamin", "Replete")))

OTB$Time_in_Days <- OTB$Time / 24

IRL water:

IRLrec<-as.data.frame(read_csv("IRL_1006.csv"))
IRLrec$Time <- as.numeric(sub('.', '', IRLrec$Time))
IRLrec$Date <- ymd(IRLrec$Date)
IRLrec$Time <- IRLrec$Time * 48
IRLrec$Time_in_Days <- IRLrec$Time / 24  # Convert Time from hours to days

IRL culture:

IRLcul <- as.data.frame(read_csv("1003_1006.csv"))
IRLcul$Time <- as.numeric(gsub("\\.", "", IRLcul$Time))  # Remove all periods
IRLcul$Date <- ymd(IRLcul$Date)
IRLcul$Time <- IRLcul$Time * 48
IRLcul$Time_in_Days <- IRLcul$Time / 24 

4.2 Combine all datasets into one with a new ‘Dataset’ column:

all_experiments <- bind_rows(
  fluor %>% mutate(Dataset = "Low Diversity OTB"),
  IRLcul %>% mutate(Dataset = "IRL culture microbiome"),
  OTB %>% mutate(Dataset = "OTB water microbiome"),
  IRLrec %>% mutate(Dataset = "IRL water microbiome")
) %>%mutate(
  Treatment = recode(Treatment,
                     "-Cobalamin" = "\u2013B12",
                     "-Biotin" = "\u2013Biotin",
                     "-Thiamine" = "\u2013Thiamine"),
  Treatment = factor(Treatment, levels = c("\u2013B12", "Replete"))
) %>%
  filter(Round %in% c(1, 2), Treatment %in% c("\u2013B12", "Replete"))

4.3 Visulaize growth of P. bahamense with recruited bacteria

# Filter out "Low Diversity OTB"
all_experiments_filtered <- all_experiments %>%
  filter(Dataset != "Low Diversity OTB") %>%
  mutate(
    Treatment = recode(Treatment,
                       "–B12" = "—B12",  # convert en dash to em dash
                       "Replete" = "Replete"),
    Treatment = factor(Treatment, levels = c("—B12", "Replete")),
    Dataset = factor(Dataset, levels = c(
      "IRL culture microbiome",
      "OTB water microbiome",
      "IRL water microbiome"
    ))
  )

growth <- ggplot(all_experiments_filtered, aes(x = Time_in_Days, y = `Chlorophyll-A (RFU)`, color = Treatment)) +
  geom_point(size = 3) +
  geom_smooth(se = FALSE) +
  scale_color_manual(values = c("—B12" = "#102E50", "Replete" = "#F5C45E")) +
  facet_grid(Round ~ Dataset) +
  labs(
    x = "Time (days)",
    y = expression(paste("Chlorophyll-", italic("a"), " fluorescence (RFU)"))
  ) +
  scale_y_continuous(labels = label_comma()) +
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 16, color = "black"),
    axis.text.x = element_text(size = 12, color = "black"),
    axis.title.y = element_text(size = 20, margin = margin(t = 10)),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 20),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 16, color = "black"),
    legend.text = element_text(size = 16, color = "black"),
    legend.position = "bottom",
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt"),
    panel.spacing = unit(.1, "cm") ) +
  coord_cartesian(ylim = c(0, 5000))

plot(growth)

Create new variable for experimemt and combine experimental data into dataframe:

fluor$Experiment <- "fluor"  
IRLcul$Experiment <- "IRLcul"  
OTB$Experiment <- "OTB"  
IRLrec$Experiment <- "IRLwat"  
all_exp <- bind_rows(fluor, IRLcul, OTB, IRLrec) %>%filter(Round %in% c(1, 2, 3,4)) %>%
  filter(Treatment %in% c("-Cobalamin", "Replete")) 

4.4 Calculate maximum RFU values for each round and medium treatment

max_bio <- all_exp %>%
  group_by(Treatment,  Bio_Rep,Experiment, Round) %>%
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE), .groups = "drop") %>%
  arrange(Treatment,  Bio_Rep)
print(max_bio)
## # A tibble: 44 × 5
##    Treatment  Bio_Rep Experiment Round Max_Chl
##    <chr>        <dbl> <chr>      <dbl>   <dbl>
##  1 -Cobalamin       1 IRLcul         1   103. 
##  2 -Cobalamin       1 IRLcul         2   136. 
##  3 -Cobalamin       1 IRLwat         1   893. 
##  4 -Cobalamin       1 IRLwat         2   904. 
##  5 -Cobalamin       1 OTB            1   341. 
##  6 -Cobalamin       1 OTB            2    87.8
##  7 -Cobalamin       2 IRLcul         1   147. 
##  8 -Cobalamin       2 IRLcul         2   173. 
##  9 -Cobalamin       2 IRLwat         1   674. 
## 10 -Cobalamin       2 IRLwat         2  1126. 
## # ℹ 34 more rows

4.4.1 Plot maximum biomass

# Define experiment order
experiment_order <- c("fluor", "IRLcul", "OTB", "IRLwat")

# Filter and prepare data for -Cobalamin only
plot_data <- max_bio %>%
  filter(Treatment == "-Cobalamin") %>%
  mutate(Experiment = factor(Experiment, levels = experiment_order))
# 
# plot_summary <- summary_stats %>%
#   filter(Treatment == "-Cobalamin") %>%
#   mutate(Experiment = factor(Experiment, levels = experiment_order))

# Build plot
max_bio_plot <- ggplot(plot_data, aes(x = Experiment, y = Max_Chl, color = Treatment)) +

  # Raw replicate points (solid)
  geom_point(size = 5, shape = 16, alpha = 1,
             position = position_jitter(width = 0.15)) +

  stat_summary(
    fun = mean,
    geom = "point",
    shape = 1,
    size = 5,
    stroke = 1.2
  ) +
  
  # Error bars: mean ± SD  (NOT SE)
  stat_summary(
    fun.data = mean_sdl,
    fun.args = list(mult = 1),   # 1 SD above/below mean
    geom = "errorbar",
    width = 0.2
  ) +

  # Legend with open circle
  scale_color_manual(
    values = c("-Cobalamin" = "#102E50"),
    labels = c("-Cobalamin" = expression("—B"[12])),
    name = "Treatment",
    guide = guide_legend(
      override.aes = list(
        shape = 21,
        fill = "white",   # open circle in legend
        stroke = 1.2,
        size = 5
      )
    )
  ) +

  # Axis labels
  labs(
    x = "Source",
    y = expression(paste("Maximum RFU"))
  ) +

  # X-axis labels
  scale_x_discrete(labels = c(
    fluor = "Low Diversity",
    IRLcul = "IRL Culture",
    OTB = "OTB Seawater",
    IRLwat = "IRL Seawater"
  )) +

  # Y-axis formatting
  scale_y_continuous(labels = label_comma()) +

  # # Facet by Round
  # facet_grid(~Round) +

  # Theme
  theme_minimal() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid = element_blank(),
    axis.text.y = element_text(size = 20, color = "black"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 14, color = "black"),
    axis.title.y = element_text(size = 20),
    axis.title.x = element_text(size = 20, margin = margin(t = 10)),
    strip.text.x = element_text(size = 20),
    strip.text = element_text(size = 20, color = "black"),
    legend.title = element_text(size = 16, color = "black"),
    legend.text = element_text(size = 16, color = "black"),
    legend.position = "bottom",
    legend.key.size = unit(1.5, "cm"),
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +

  # Set visible y-axis range safely
  coord_cartesian(ylim = c(50, 1500))

# Print plot
max_bio_plot

4.4.2 Run an ANOVA followed by a Tukey HSD test for the experiment factor:

max_B12 <- all_exp %>%
  group_by(Experiment ,Treatment, Round, Bio_Rep) %>% 
  summarise(Max_Chl = max(`Chlorophyll-A (RFU)`, na.rm = TRUE)) %>%
  ungroup()  %>% filter(Treatment %in% c("-Cobalamin"))


aov_biomass <- aov(Max_Chl ~ Experiment, data =max_B12)
summary(aov_biomass)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## Experiment   3 2481763  827254   38.95 4.45e-08 ***
## Residuals   18  382327   21240                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov_biomass)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Max_Chl ~ Experiment, data = max_B12)
## 
## $Experiment
##                     diff       lwr       upr     p adj
## IRLcul-fluor  -148.55167 -414.4357  117.3324 0.4145013
## IRLwat-fluor   679.19000  413.3060  945.0740 0.0000057
## OTB-fluor      -23.03833 -288.9224  242.8457 0.9946502
## IRLwat-IRLcul  827.74167  589.9278 1065.5556 0.0000001
## OTB-IRLcul     125.51333 -112.3006  363.3272 0.4626003
## OTB-IRLwat    -702.22833 -940.0422 -464.4144 0.0000007

4.5 Specific growth rates using all_easylinear

4.5.1 Remove outliers from IRL culture data:

outliers are considered values the drop below T0 RFU values

# Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(IRLcul, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
IRLcul_merged <- merge(IRLcul, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
IRLcul_clean <- subset(
  IRLcul_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

4.5.2 Calculate specific growth rates for IRL culture:

IRL culture trial 1 Replete:

IRLcul_mumax <- IRLcul_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)
IRLcul1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRLcul_mumax, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLcul1 , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLcul1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.03861  0.16034 -0.15112  0.16012 -0.18606 -0.07377  0.12911 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.7592193  0.3438351   5.116  0.00372 ** 
## x           0.0075785  0.0006407  11.828  7.6e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1627 on 5 degrees of freedom
## Multiple R-squared:  0.9655, Adjusted R-squared:  0.9586 
## F-statistic: 139.9 on 1 and 5 DF,  p-value: 7.603e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.05977  0.09664  0.02457  0.12515 -0.33697  0.07566  0.07473 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.793494   0.469100   1.692 0.151524    
## x           0.007143   0.000691  10.336 0.000146 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1755 on 5 degrees of freedom
## Multiple R-squared:  0.9553, Adjusted R-squared:  0.9463 
## F-statistic: 106.8 on 1 and 5 DF,  p-value: 0.0001459
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.03760 -0.17761  0.53725 -0.39227  0.01679  0.10783 -0.05440 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.488243   0.895579   1.662  0.15745   
## x           0.005008   0.001233   4.062  0.00971 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3132 on 5 degrees of freedom
## Multiple R-squared:  0.7674, Adjusted R-squared:  0.7209 
## F-statistic:  16.5 on 1 and 5 DF,  p-value: 0.009713
# Extracts only the estimated regression coefficients
coef(IRLcul1)
##                y0    y0_lm       mumax      lag
## 1:Replete:1 74.22 5.807901 0.007578469 336.1912
## 1:Replete:2 88.61 2.211109 0.007142610 516.7230
## 1:Replete:3 71.52 4.429309 0.005007900 555.4690
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLcul1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9654963      0.9552895      0.7674161

IRL culture trial 2 Replete:

IRLcul_mumax2 <- IRLcul_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLcul2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRLcul_mumax2, h =7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLcul2 , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLcul2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.23499  0.44031 -0.16656 -0.01381 -0.09558  0.10722 -0.03659 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.7712453  0.4712913   5.880 0.002020 ** 
## x           0.0069299  0.0009628   7.198 0.000806 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2445 on 5 degrees of freedom
## Multiple R-squared:  0.912,  Adjusted R-squared:  0.8944 
## F-statistic: 51.81 on 1 and 5 DF,  p-value: 0.0008063
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.15234 -0.55703  0.04284  0.55888  0.12964 -0.23617 -0.09050 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.626180   0.806476   2.016  0.09983 . 
## x           0.008378   0.001503   5.575  0.00256 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3817 on 5 degrees of freedom
## Multiple R-squared:  0.8614, Adjusted R-squared:  0.8337 
## F-statistic: 31.08 on 1 and 5 DF,  p-value: 0.002558
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.31177 -0.51835 -0.13281  0.31474  0.18935 -0.07054 -0.09416 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 2.906698   0.628814   4.623  0.00572 **
## x           0.006738   0.001285   5.245  0.00334 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3263 on 5 degrees of freedom
## Multiple R-squared:  0.8462, Adjusted R-squared:  0.8155 
## F-statistic: 27.51 on 1 and 5 DF,  p-value: 0.003339
# Extracts only the estimated regression coefficients
coef(IRLcul2)
##                y0     y0_lm       mumax      lag
## 2:Replete:1 73.64 15.978520 0.006929859 220.4869
## 2:Replete:2 81.69  5.084417 0.008378133 331.4284
## 2:Replete:3 71.42 18.296276 0.006738242 202.1121
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLcul2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9119822      0.8614255      0.8462232

4.5.3 Remove outliers from OTB water data:

# Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(OTB, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
OTBwat_merged <- merge(OTB, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
OTBwat_clean <- subset(
  OTBwat_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

4.5.4 Calculate specific growth rates for OTB water:

OTB water trial 1 Replete:

OTB_mumax1 <- OTBwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(1),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


OTBwat1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, OTB_mumax1, h = 7, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(OTBwat1 , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(OTBwat1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## -0.15859  0.22093 -0.08966 -0.03358  0.16311 -0.01998 -0.08224 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.332309   0.210366   20.59 5.00e-06 ***
## x           0.007191   0.000602   11.95 7.25e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1529 on 5 degrees of freedom
## Multiple R-squared:  0.9661, Adjusted R-squared:  0.9594 
## F-statistic: 142.7 on 1 and 5 DF,  p-value: 7.251e-05
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.33247 -0.31312 -0.07991 -0.10915  0.06674  0.08438  0.01858 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.6685656  0.3405662  10.772 0.000120 ***
## x           0.0084842  0.0008604   9.861 0.000183 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2185 on 5 degrees of freedom
## Multiple R-squared:  0.9511, Adjusted R-squared:  0.9413 
## F-statistic: 97.23 on 1 and 5 DF,  p-value: 0.0001828
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  0.05109  0.05233 -0.25963  0.01366  0.20820  0.01297 -0.07862 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.3310777  0.2440626   17.75 1.04e-05 ***
## x           0.0068159  0.0006166   11.05 0.000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1566 on 5 degrees of freedom
## Multiple R-squared:  0.9607, Adjusted R-squared:  0.9528 
## F-statistic: 122.2 on 1 and 5 DF,  p-value: 0.0001055
# Extracts only the estimated regression coefficients
coef(OTBwat1)
##                 y0    y0_lm       mumax      lag
## 1:Replete:1 211.13 76.11986 0.007190613 141.8745
## 1:Replete:2 101.61 39.19564 0.008484247 112.2759
## 1:Replete:3 328.76 76.02618 0.006815913 214.8282
# Returns the model’s R² (proportion of variance explained)
rsquared(OTBwat1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9661413      0.9510923      0.9606888

OTB water trial 2 Replete:

OTB_mumax2 <- OTBwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>% 
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)


OTBwat2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, OTB_mumax2, h = 5, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(OTBwat2  , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(OTBwat2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.04205 -0.14349  0.12612  0.01002 -0.03470 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.0096618  0.0889572  56.315 1.23e-05 ***
## x           0.0041960  0.0007566   5.546   0.0116 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1148 on 3 degrees of freedom
## Multiple R-squared:  0.9111, Adjusted R-squared:  0.8815 
## F-statistic: 30.76 on 1 and 3 DF,  p-value: 0.01156
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.12058  0.06702 -0.35136  0.01932  0.14443 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.724234   0.180952  26.108 0.000123 ***
## x           0.007660   0.001539   4.977 0.015586 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2336 on 3 degrees of freedom
## Multiple R-squared:  0.892,  Adjusted R-squared:  0.856 
## F-statistic: 24.77 on 1 and 3 DF,  p-value: 0.01559
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##       1       2       3       4       5 
##  0.1065  0.1344 -0.4463  0.0634  0.1420 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.409273   0.224859  19.609  0.00029 ***
## x           0.007442   0.001912   3.891  0.03009 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2903 on 3 degrees of freedom
## Multiple R-squared:  0.8347, Adjusted R-squared:  0.7795 
## F-statistic: 15.14 on 1 and 3 DF,  p-value: 0.03009
# Extracts only the estimated regression coefficients
coef(OTBwat2)
##                 y0     y0_lm       mumax      lag
## 2:Replete:1 156.29 149.85405 0.004195997 10.02180
## 2:Replete:2 127.08 112.64419 0.007660107 15.74165
## 2:Replete:3  91.45  82.20965 0.007442250 14.31284
# Returns the model’s R² (proportion of variance explained)
rsquared(OTBwat2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.9111292      0.8919806      0.8346501

4.5.5 Remove outliers from IRL water data:

#Compute T0 baseline per replicate, treatment, and round
baseline <- aggregate(
  `Chlorophyll-A (RFU)` ~ Bio_Rep + Treatment + Round,
  data = subset(IRLrec, Time == 0),
  FUN = mean,
  na.rm = TRUE
)
names(baseline)[names(baseline) == "Chlorophyll-A (RFU)"] <- "Chl_T0"

# Merge baseline back into the full dataset
IRLwat_merged <- merge(IRLrec, baseline, by = c("Bio_Rep", "Treatment", "Round"), all.x = TRUE)

# Keep all T0 points, and remove any post-T0 points where replicate < its T0
IRLwat_clean <- subset(
  IRLwat_merged,
  Time == 0 | `Chlorophyll-A (RFU)` >= Chl_T0
)

4.5.6 Calculate specific growth rates for IRL water:

IRL water trial 1 Replete:

IRL_mumax1 <- IRLwat_clean %>%
  filter(Treatment %in% c("Replete"), Round %in% c(1)) %>%
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  droplevels()

IRLwat1 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax1, h = 5, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat1 , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLwat1)
## $`1:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.07587  0.09142 -0.13804 -0.22946  0.35195 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.005693   0.268748  18.626 0.000338 ***
## x           0.008304   0.001006   8.258 0.003718 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2644 on 3 degrees of freedom
## Multiple R-squared:  0.9579, Adjusted R-squared:  0.9438 
## F-statistic:  68.2 on 1 and 3 DF,  p-value: 0.003718
## 
## 
## $`1:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.01569 -0.18763  0.22688 -0.29711  0.27355 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.6061505  0.2231989   25.12 0.000138 ***
## x           0.0064232  0.0008246    7.79 0.004403 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2887 on 3 degrees of freedom
## Multiple R-squared:  0.9529, Adjusted R-squared:  0.9372 
## F-statistic: 60.68 on 1 and 3 DF,  p-value: 0.004403
## 
## 
## $`1:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.32328  0.13706  0.35205 -0.08458 -0.08125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.230102   0.228469  18.515 0.000344 ***
## x           0.013307   0.001943   6.848 0.006374 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.295 on 3 degrees of freedom
## Multiple R-squared:  0.9399, Adjusted R-squared:  0.9198 
## F-statistic: 46.89 on 1 and 3 DF,  p-value: 0.006374
# Extracts only the estimated regression coefficients
coef(IRLwat1)
##                 y0     y0_lm       mumax        lag
## 1:Replete:1 108.14 149.26049 0.008304417 -38.806618
## 1:Replete:2 267.86 272.09478 0.006423229  -2.442075
## 1:Replete:3 178.44  68.72424 0.013306580  71.705157
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLwat1)
## 1:Replete:1.r2 1:Replete:2.r2 1:Replete:3.r2 
##      0.9578639      0.9528881      0.9398718

IRL water trial 2 Replete:

IRL_mumax2 <- IRLwat_clean %>%
  filter(Treatment %in% c("Replete"),
         Round %in% c(2),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat2 <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax2, h =6, quota = 1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat2 , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLwat2)
## $`2:Replete:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
##  0.27656 -0.39465 -0.12399  0.08819  0.39137 -0.23748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.126211   0.317395  13.000 0.000202 ***
## x           0.009866   0.001698   5.811 0.004365 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3409 on 4 degrees of freedom
## Multiple R-squared:  0.8941, Adjusted R-squared:  0.8676 
## F-statistic: 33.76 on 1 and 4 DF,  p-value: 0.004365
## 
## 
## $`2:Replete:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
##  0.03290  0.04031 -0.06486  0.15651 -0.44416  0.27930 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.883134   0.200405  24.366 1.68e-05 ***
## x           0.006915   0.001379   5.014  0.00741 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2769 on 4 degrees of freedom
## Multiple R-squared:  0.8628, Adjusted R-squared:  0.8284 
## F-statistic: 25.14 on 1 and 4 DF,  p-value: 0.007414
## 
## 
## $`2:Replete:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5        6 
## -0.10408  0.24742  0.10827 -0.59427  0.39445 -0.05179 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.516376   0.359111  12.577  0.00023 ***
## x           0.008128   0.001921   4.231  0.01336 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3857 on 4 degrees of freedom
## Multiple R-squared:  0.8174, Adjusted R-squared:  0.7717 
## F-statistic:  17.9 on 1 and 4 DF,  p-value: 0.01336
# Extracts only the estimated regression coefficients
coef(IRLwat2)
##                 y0     y0_lm       mumax       lag
## 2:Replete:1 102.18  61.94276 0.009865993 50.732380
## 2:Replete:2 136.46 132.04390 0.006914858  4.757449
## 2:Replete:3 105.76  91.50335 0.008128438 17.813607
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLwat2)
## 2:Replete:1.r2 2:Replete:2.r2 2:Replete:3.r2 
##      0.8940791      0.8627528      0.8173779

IRL water trial 1 -B12:

IRL_mumax1_b <- IRLwat_clean %>%
  filter(Treatment %in% c("-Cobalamin"),
         Round %in% c(1),
         Time >= 0) %>%  # <-- exclude T0
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat1_b <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax1_b, h = 4, quota =1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat1_b , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLwat1_b)
## $`1:-Cobalamin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
## -0.13839  0.02356  0.36806 -0.25322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.629046   0.277043  16.709  0.00356 **
## x           0.011942   0.003085   3.871  0.06073 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3311 on 2 degrees of freedom
## Multiple R-squared:  0.8822, Adjusted R-squared:  0.8233 
## F-statistic: 14.98 on 1 and 2 DF,  p-value: 0.06073
## 
## 
## $`1:-Cobalamin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.09071 -0.12902 -0.06650  0.10480 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.1181197  0.1239394  41.295 0.000586 ***
## x           0.0053753  0.0007698   6.982 0.019900 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1419 on 2 degrees of freedom
## Multiple R-squared:  0.9606, Adjusted R-squared:  0.9409 
## F-statistic: 48.76 on 1 and 2 DF,  p-value: 0.0199
## 
## 
## $`1:-Cobalamin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4 
##  0.06389 -0.33062  0.46957 -0.20284 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 3.332327   0.711542   4.683   0.0427 *
## x           0.015666   0.004035   3.883   0.0604 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.433 on 2 degrees of freedom
## Multiple R-squared:  0.8829, Adjusted R-squared:  0.8243 
## F-statistic: 15.08 on 1 and 2 DF,  p-value: 0.06038
# Extracts only the estimated regression coefficients
coef(IRLwat1_b)
##                    y0     y0_lm       mumax       lag
## 1:-Cobalamin:1  89.18 102.41631 0.011941636 -11.58880
## 1:-Cobalamin:2 182.88 167.02102 0.005375287  16.87548
## 1:-Cobalamin:3 112.12  28.00342 0.015665946  88.55150
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLwat1_b)
## 1:-Cobalamin:1.r2 1:-Cobalamin:2.r2 1:-Cobalamin:3.r2 
##         0.8822319         0.9605951         0.8828865

IRL water trial 2 -B12:

IRL_mumax2_b <-  IRLwat_clean %>%
  filter(Treatment %in% c("-Cobalamin"),
         Round %in% c(2),
         Time >= 0) %>%  
  rename(Chl = `Chlorophyll-A (RFU)`) %>%
  arrange(Time, Bio_Rep)

IRLwat2_b <- all_easylinear(Chl ~ Time|Round+Treatment+Bio_Rep, IRL_mumax2_b, h =5, quota =1)

par(mfrow = c(1, 3))
par(mar = c(2, 2, 2, 1))
plot(IRLwat2_b , log = "y")

# Shows the full regression output including coefficients, standard errors, t-tests, p-values, and model fit statistics
summary(IRLwat2_b)
## $`2:-Cobalamin:1`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.08915 -0.35337  0.17809  0.34733 -0.26120 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.371518   0.460476   9.493  0.00248 **
## x           0.008698   0.002261   3.847  0.03101 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3432 on 3 degrees of freedom
## Multiple R-squared:  0.8314, Adjusted R-squared:  0.7752 
## F-statistic:  14.8 on 1 and 3 DF,  p-value: 0.03101
## 
## 
## $`2:-Cobalamin:2`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
## -0.04172  0.31220 -0.17814 -0.37176  0.27941 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.022854   0.290615  17.284 0.000422 ***
## x           0.007184   0.001842   3.899 0.029934 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3402 on 3 degrees of freedom
## Multiple R-squared:  0.8352, Adjusted R-squared:  0.7803 
## F-statistic:  15.2 on 1 and 3 DF,  p-value: 0.02993
## 
## 
## $`2:-Cobalamin:3`
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##        1        2        3        4        5 
##  0.23747  0.05598 -0.51117 -0.09547  0.31319 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 4.206086   0.396161  10.617  0.00179 **
## x           0.010498   0.002488   4.219  0.02434 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3777 on 3 degrees of freedom
## Multiple R-squared:  0.8558, Adjusted R-squared:  0.8077 
## F-statistic:  17.8 on 1 and 3 DF,  p-value: 0.02434
# Extracts only the estimated regression coefficients
coef(IRLwat2_b)
##                    y0     y0_lm       mumax       lag
## 2:-Cobalamin:1 107.43  79.16373 0.008697853 35.103053
## 2:-Cobalamin:2 145.64 151.84401 0.007183847 -5.806909
## 2:-Cobalamin:3 130.09  67.09340 0.010498480 63.070164
# Returns the model’s R² (proportion of variance explained)
rsquared(IRLwat2_b)
## 2:-Cobalamin:1.r2 2:-Cobalamin:2.r2 2:-Cobalamin:3.r2 
##         0.8314296         0.8352092         0.8557594

4.5.7 Extract and combine all specific growth rates:

# Helper function to handle matrix or vector coefficients
make_df <- function(model, round_num, experiment_name) {
  coefs <- coef(model)
  
  if (is.matrix(coefs)) {
    df <- as.data.frame(coefs) %>%
      tibble::rownames_to_column("Group") %>%
      mutate(
        Round = round_num,
        Experiment = experiment_name,
        r2 = rsquared(model)
      )
  } else {
    df <- data.frame(
      Group = names(coefs),
      Estimate = as.numeric(coefs),
      Round = round_num,
      Experiment = experiment_name,
      r2 = rep(rsquared(model), length(coefs))
    )
  }
  
  return(df)
}

# Create data frames for all models with experiment names
df1 <- make_df(IRLcul1, 1, "IRL Culture")
df2 <- make_df(IRLcul2, 2, "IRL Culture")
df3 <- make_df(OTBwat1, 1, "OTB Water")
df4 <- make_df(OTBwat2, 2, "OTB Water")
df5 <- make_df(IRLwat1, 1, "IRL Water")
df6 <- make_df(IRLwat2, 2, "IRL Water")
df7 <- make_df(IRLwat1_b, 1, "IRL Water B")
df8 <- make_df(IRLwat2_b, 2, "IRL Water B")

# Combine all results
all_results <- bind_rows(df1, df2, df3, df4, df5, df6, df7, df8)

# Optional: reorder columns
all_results <- all_results %>%
  select(Experiment, Round, Group, y0, mumax, r2)

4.5.8 Prepare IRL seawater data for T-test for both trials in -B12 media:

df5 <- as.data.frame(coef(IRLwat1))
df5$Round <- 1
df5$Experiment <- "IRL Water"

df6 <- as.data.frame(coef(IRLwat2))
df6$Round <- 2
df6$Experiment <- "IRL Water"

df7 <- as.data.frame(coef(IRLwat1_b))
df7$Round <- 1
df7$Experiment <- "IRL Water -B12"

df8 <- as.data.frame(coef(IRLwat2_b))
df8$Round <- 2
df8$Experiment <- "IRL Water -B12"

# Combine
IRLwater <- bind_rows(df5, df6, df7, df8)

T-test between specific growthrates both trials for IRL water in -B12:

t.test(mumax ~ Experiment, data = IRLwater, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  mumax by Experiment
## t = -0.59074, df = 10, p-value = 0.5678
## alternative hypothesis: true difference in means between group IRL Water and group IRL Water -B12 is not equal to 0
## 95 percent confidence interval:
##  -0.005105451  0.002965605
## sample estimates:
##      mean in group IRL Water mean in group IRL Water -B12 
##                  0.008823919                  0.009893842

4.5.9 Plot IRL seawater specific growth rates:

# Extract the columns needed
mumax_results <- all_results %>%
  select(Experiment, Round, Group, mumax) %>%
  mutate(ExperimentRound = paste(Experiment, Round, sep = "_"))

mumax_results <- mumax_results %>%
  mutate(
    Treatment = factor(
      Experiment,
      levels = c("IRL Water B", "IRL Water")  # order: B12 first, Replete second
    )
  )

# Filter IRL Water only
mumax_IRLwater <- mumax_results %>%
  filter(grepl("IRL Water", Experiment))


IRL_mu <- ggplot(mumax_IRLwater, aes(x = Treatment, y = mumax, color = Treatment)) +
  # Raw replicate points
  geom_point(size = 5, shape = 16, position = position_jitter(width = 0.1)) +
  
  # Mean per Experiment+Round (open circle)
  stat_summary(fun = mean, geom = "point", shape = 1, size = 5, stroke = 1.2) +
  
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  
  # facet_grid(~Round) +
  scale_color_manual(values = c(
    "IRL Water B" = "#102E50",
    "IRL Water" = "#F5C45E",
    "-Thiamine" = "deeppink4",
    "-Biotin" = "#667028"
  )) +
  scale_x_discrete(labels = c(
  "IRL Water B"  = expression("\u2013B"[12]), 
  "IRL Water"   = "Replete",
  "-Thiamine"   = expression("\u2013Thiamine"),
  "-Biotin"     = expression("\u2013Biotin")    
))+
  theme_minimal()  +

  theme(
    legend.position = "none",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "black"),
    axis.title = element_text(size = 20),
    axis.text = element_text(size = 20, color = "black"),
    axis.text.x = element_text( size = 18, color = "black", angle = 45, hjust = 1),
    strip.text = element_text( size = 20, color = "black"),
    legend.title = element_blank(),
    legend.text = element_blank(),
    legend.key = element_blank(),
    axis.ticks = element_line(color = "black", linewidth = 0.5),
    axis.ticks.length = unit(4, "pt")
  ) +
  labs(
    x = "Media Treatment",
    y = expression(paste("(", mu, ") RFU" * "\u00B7" * "day"^-1, " "))
  ) +
  guides(color = "none", fill = "none", shape = "none") +
  
  # Set visible range safely
  coord_cartesian(ylim = c(0, 0.5)) 
  
IRL_mu

4.6 Figure 5

# Ensure BOTH subplots have no legend inside patchwork
max_bio__2 <- max_bio_plot + theme(legend.position = "none")
mu_2      <- IRL_mu + theme(legend.position = "none")

bo <- 
  growth / 
  (max_bio_plot  +
   IRL_mu ) + 
  plot_layout(
    guides = "collect",
    widths = c(4, 2)     
  ) &
  theme(
    legend.position = "bottom",
    panel.spacing = unit(.1, "lines"),
    plot.margin = margin(15, 15, 15, 15),
    axis.title.y = element_text(margin = margin(r = 12)),
    axis.title.x = element_text(margin = margin(t = 8)),
    legend.margin = margin(t = 4, b = 4),
    legend.box.margin = margin(5, 5, 5, 5),
    axis.text.x = element_text(size = 16)
  )

bo

Fig. 5 Growth responses of P. bahamense (FWC1006, OTB) with bacterial communities recruited from FWC1003 (IRL) P. bahamense culture, OTB seawater, and IRL seawater when grown in L1/4 media missing vitamin B12 and replete L1/4 media A In vivo Chlorophyll a fluorescence (Chl-a) was used to track the growth of P. bahamense (FWC1006, OTB) with bacterial communities recruited from three sources in L1/4 media missing B12 (yellow) and replete L1/4 media (blue) for two experimental rounds. Replicates growing in replete media and replicates with bacteria recruited from IRL seawater were serially transferred from experimental round one to round two. The biomass in replicates with bacteria from the FWC1003 P. bahamense culture and from OTB seawater grown in media missing B12 were too low to support serial transfers; the replicates in these treatments were restarted from cultures with the same bacterial communities in replete media for round two. B Maximum Chl-a values in each replicate grown in media missing B12 by bacterial community source. The maximum Chl-a reached in media missing B12 with the initial low-diversity bacterial community in the FWC1006 P. bahamense culture (from figure 3) is included for reference. One-way ANOVA results indicate that bacterial community source had a significant effect on maximum Chl-a when vitamin B12 was excluded from media (p = 0.0001, F = 50.27). Post-hoc multiple comparisons (Tukey HSD) showed that the cultures with bacteria from IRL seawater had significantly higher maximum Chl-a compared to cultures with all other bacterial communities (p-adj = 0.0000 for each comparison). The maximum Chl-a was not significantly different between cultures with the other three bacterial communities. Results from key comparisons are shown on the plot with the following significance codes: p-adj < 0.001; ’*’. C** Specific growth rates (µ; RFU day–1) for each replicate with the bacterial community recruited from IRL seawater by media formulation. Means are plotted as black open circles, and error bars represent one standard deviation of the mean. The growth rates for FWC 1006 P. bahamense with bacteria from IRL seawater were not significantly different in media with or without B12 (two-sample t-test; p = 0.5678, t(10) = –0.5907).

5 Characterizing P. bahamense recruited bacterial communities

  1. Full-length 16S rRNA was sequenced using Oxford Nanopore MinIon v10.4 flow cells (Oxford Nanopore Technology 2008)
  2. Sequence data was processed and classified using Kraken v2.1.2 (Lu and Salzberg 2020)

5.1 Import data:

Recruit_16s <- read.csv("/Users/lydiaruggles/Desktop/KRAKEN-16s-counts-genus.csv")

5.2 Convert to long data:

# Identify sample columns
sample_cols <- names(Recruit_16s)[2:(which(names(Recruit_16s) == "total") - 1)]

# Pivot to long data
Recruit_16s_long <- Recruit_16s %>%
  pivot_longer(
    cols = all_of(sample_cols),
    names_to = "sample",
    values_to = "count"
  )

# Separate samples into new columns
Recruit_16s_long <- Recruit_16s_long %>%
  separate(sample, into = c("experiment","treatment","filter_size","bio_rep"),
           sep = "_", remove = FALSE)

# One genus in one sample with the total count of that genus in that sample
grouped_data1 <- Recruit_16s_long %>%
  mutate(count = as.numeric(count)) %>%
  group_by(treatment, experiment, filter_size, bio_rep, genus) %>%
  summarise(total_count = sum(count), .groups = "drop"
)

5.3 Prepare data as wide matrix:

# Pivot wider
wide_data1 <- grouped_data1 %>%
  pivot_wider(names_from = genus, values_from = total_count, values_fill = 0)

# Create rownames and remove metadata columns for distance matrix
rownames(wide_data1) <- paste(wide_data1$treatment, wide_data1$experiment, wide_data1$filter_size, wide_data1$bio_rep, sep = "_")

wide_matrix1 <- wide_data1[, !(names(wide_data1) %in% c("treatment", "experiment", "filter_size", "bio_rep"))]
wide_matrix1 <- as.data.frame(lapply(wide_matrix1, as.numeric))

wide_matrix1[is.na(wide_matrix1)] <- 0

5.4 Compute Bray–Curtis Distances and Perform PCoA

# Bray-Curtis + PCoA
bray_dist1 <- vegdist(wide_matrix1, method = "bray")
pcoa_result1 <- pcoa(bray_dist1)

# Build plot data
pcoa_data1 <- as.data.frame(pcoa_result1$vectors)
pcoa_data1$treatment <- wide_data1$treatment
pcoa_data1$experiment <- wide_data1$experiment
pcoa_data1$filter_size <- as.factor(wide_data1$filter_size)

5.5 Calculate Permanova:

#ensure wide matrix is numeric
wide_matrix1 <- as.data.frame(lapply(wide_matrix1, as.numeric))
wide_matrix1[is.na(wide_matrix1)] <- 0

#Build the Bray–Curtis distance matrix
bray_dist <- vegdist(wide_matrix1, method = "bray")

#Bring metadata back in (from wide_data1)
metadata <- wide_data1[, c("treatment", "experiment", "filter_size", "bio_rep")]


permanova_result <- adonis2(
  bray_dist ~ experiment,
  data = metadata,
  permutations = 999
)

print(permanova_result)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = bray_dist ~ experiment, data = metadata, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## Model     4   6.2820 0.67542 16.127  0.001 ***
## Residual 31   3.0189 0.32458                  
## Total    35   9.3009 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.6 Figure 6 PCoA plot

sw_pcoa <- ggplot(
  pcoa_data1,
  aes(
    x = Axis.1,
    y = Axis.2,
    color = experiment,
    shape = treatment,
    alpha = filter_size   # keep as factor
  )
) +
  geom_point(size = 5) +
  labs(
    x = paste0("PCoA1 (", round(pcoa_result1$values$Relative_eig[1] * 100, 2), "%)"),
    y = paste0("PCoA2 (", round(pcoa_result1$values$Relative_eig[2] * 100, 2), "%)"),
    color = "Experiment",
    shape = "Treatment",
    alpha = "Filter size"
  ) +
  stat_ellipse(aes(group = experiment)) +
  scale_color_manual(values = pp_palette)+

  # Map discrete alpha manually
  scale_alpha_manual(values = c("02" = 0.4, "10" = 1)) +

  theme_minimal() +
  theme(
    legend.key.size = unit(1.5, "cm"),
    legend.title = element_text(size = 16),
    legend.text = element_text(size =12),
    panel.border = element_rect(color = 'black', fill = NA),
    axis.title.x = element_text(size = 16, color = "black"),
    axis.title.y = element_text(size =16, color = "black"),
    strip.text = element_text(size =12, color = "black"),
    axis.text.y = element_text(size = 16, color = "black"),
    axis.text.x = element_text(size = 16, color = "black")
  ) +
  geom_hline(yintercept = 0, linetype = "dotted", linewidth = 0.5) +
  geom_vline(xintercept = 0, linetype = "dotted", linewidth = 0.5) +
  theme(panel.grid = element_blank())

plot(sw_pcoa)

5.7 Prepare data for relative abundance plot:

#Standardize experiment names
grouped_data <- grouped_data1 %>%
  mutate(
    experiment = case_when(
      grepl("1003.1006", experiment, ignore.case = TRUE) ~ "1003.1006",
      grepl("OTB", experiment, ignore.case = TRUE) ~ "OTB.1006",
      grepl("IRL", experiment, ignore.case = TRUE) ~ "IRL.1006",
      TRUE ~ experiment
    ),
    experiment = factor(experiment,
                        levels = c("1003.1006", "OTB_seawater", "OTB.1006", "IRL_seawater", "IRL.1006"),
                        labels = c("Microbiome from IRL culture",
                                   "OTB Seawater",
                                   "Microbiome from OTB seawater",
                                   "IRL seawater",
                                   "Microbiome from IRL seawater"))
  )

#Pivot to long format
sample_cols <- names(Recruit_16s)[2:(which(names(Recruit_16s) == "total") - 1)]

Recruit_16s_long <- Recruit_16s %>%
  pivot_longer(
    cols = all_of(sample_cols),
    names_to = "sample",
    values_to = "count"
  ) %>%
  separate(sample, into = c("experiment","treatment","filter_size","bio_rep"),
           sep = "_", remove = FALSE) %>%
  mutate(count = as.numeric(count))  # ensure numeric

# ️ Summarize counts per sample & genus
grouped_data <- Recruit_16s_long %>%
  group_by(experiment, treatment, filter_size, bio_rep, order) %>%
  summarise(total_count = sum(count), .groups = "drop") %>%
  group_by(experiment, treatment, filter_size, bio_rep) %>%
  mutate(RA = (total_count / sum(total_count)) * 100) %>%
  ungroup() %>%
  filter(treatment %in% c("noB12", "replete", "seawater"))

# ⃣Collapse low abundance genera (<1%) into "z<1% abund."
grouped_data$order <- as.character(grouped_data$order)
grouped_data$order [grouped_data$RA < 1] <- "z<1% abund."
grouped_data$order <- factor(grouped_data$order )  # back to factor

#  Optional: fix experiment & filter labels for plotting
grouped_data <- grouped_data %>%
  mutate(
    experiment = factor(experiment,
                        levels = c("X1003.1006","OTB","OTB.1006","IRL","IRL.1006"),
                        labels = c("Recruited IRL culture",
                                   "OTB Seawater",
                                   "Recruited OTB seawater",
                                   "IRL Seawater",
                                   "Recruited IRL seawater")),
    Treatment = recode(treatment, "-B12" = "\u2014B12"),
    filter_size = recode(filter_size, "10" = "10.0um", "02" = "0.2um"),
    filter_size = factor(filter_size, levels = c("10.0um", "0.2um")),  # <-- top first
    x_label = paste(Treatment, "Rep", bio_rep, sep = " ")
  )

#  Define color palette (ensure enough colors for all genera)
unique_genera <- levels(grouped_data$order )
n_colors <- length(unique_genera)
palette_colors <- brewer.pal(min(12, n_colors), "Set3")  # adjust if >12 genera
if(n_colors > 12) {
  palette_colors <- colorRampPalette(palette_colors)(n_colors)
}
names(palette_colors) <- unique_genera

5.8 Figure 7 Relative abundance plot

n_orders <- length(unique(grouped_data$order))

# expand the 5-color palette into n_orders colors
order_palette <- colorRampPalette(pp_palette)(n_orders)

taxa_bar_plot <- ggplot(grouped_data,
                        aes(x = x_label, y = RA, fill = order)) +
  geom_bar(stat = "identity", position = "stack") +
  facet_grid(filter_size ~ experiment, scales = "free_x") +
  theme_classic() +
  theme(
    legend.position = "right",
    legend.title = element_blank(),
    axis.text.x = element_text(size = 18, angle = 45, hjust = 1, color = "black"),
    axis.text.y = element_text(size = 18, color = "black"),
    axis.title.y = element_text(size = 18),
    strip.text = element_text(size = 12),
    legend.text = element_text(size = 18),
    strip.position = "top"
  ) +
  scale_fill_manual(values = order_palette) +
  guides(fill = guide_legend(ncol = 2)) +
  ylab("Relative Abundance (%)") +
  labs(x = NULL)

taxa_bar_plot